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CHAPTER  1 
INTRODUCTION 

When  a  turbulent  fluid  interacts  with  a  sinusoidal  solid  wave 
spatial  variations  of  the  pressure  and  the  shear  stress  occur  at  the 
surface.  If  the  wave  is  of  small  enough  amplitude  a  linear  response 
can  be  expected  in  that  the  hydrodynamic  quantities  can  be  described 
by  single  harmonics  with  characteristic  phases  and  amplitudes.  The 
analysis  for  this  case  is,  therefore,  greatly  simplified.  The  principal 
theoretical  problem  is  the  determination  of  the  Reynolds  stresses  close 
to  the  surface.  In  previous  studies  by  Cook  [15]  and  Thorsness  [70]  it 
was  shown  that  the  phase  angle  associated  with  the  spatial  variation  of 
the  surface  shear  stress  provides  a  particularly  sensitive  test  of  the 
theory  used. 

The  primary  purpose  of  this  thesis  is  to  obtain  surface  shear  stress 
measurements  for  turbulent  flow  over  small  amplitude  solid  waves,  over  a 
wide  range  of  conditions  and  to  compare  these  measurements  in  a  meaningful 
way  with  predictions  derived  from  various  turbulence  models.  These 
measurements  provide  a  particularly  sensitive  test  of  current  turbulence 
models  in  the  neighborhood  of  a  solid  boundary.  The  results  find  appli¬ 
cation  in  understanding  the  generation  of  water  waves  at  a  gas  liquid 
interface. 

Various  attempts  have  been  made  to  measure  surface  shear  stress  and 
pressure  profiles  over  a  solid  wavy  surface.  Motzfield  [49],  Larras  and 
Claria  [36]  and  Zagustin  et  al.  [76]  measured  pressure  profiles.  Their 
measurements  indicate  a  linear  response  for  amplitude,  a^,  to  wavelength,  X, 
ratios  of2a^/X  ^  0.05.  Kendall  [34],  who  was  concerned  mainly  with  moving 
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waves,  made  a  limited  number  of  pressure  and  shear  stress  measurements 
at  zero  wave  speed.  A  wave  of  X  =  4  inches  and  2aJ\=  0.062  was  used  in 
this  study.  Sigal  [65],  using  two  geometrically  similar  waves  with 

2  a  /X  =  0.055  and  wavelengths  of  6  and  12  inches,  measured  both  surface 

d 

shear  stress  and  pressure  profiles.  Hsu  and  Kennedy  [26]  carried  out  a 

similar  set  of  experiments  to  Sigal  using  waves  of  2  a d/X  -  0.022  and 

0.044  on  the  wall  of  a  pipe.  In  this  laboratory  Zilker  [77]  and  Cook  [1--] 

used  electrochemical  techniques  to  determine  shear  stress  profiles  over 

waves  with  2  a J\  =  0.0312,  0.05,  0.125  and  0.2  and  a  wavelength  of 
d 

X  =  2  inches.  The  major  limitation  of  the  above  measurements,  with  the 
exception  of  the  shear  stress  measurements  of  Zilker  and  Cook,  is  that 
they  were  not  carried  out  over  a  large  enough  range  of  flow  conditions 
to  be  suitable  for  testing  solutions  of  the  momentum  equations.  Another 
limitation  is  that  the  height  of  the  waves  was  such  that  the  appearance 
of  higher  order  harmonics  in  the  shear  stress  profiles  made  comparison 
of  these  measurements  with  linear  theory  difficult. 

Thorsness  [70]  and  Morrisroe  [48]  obtained  a  set  of  shear  stress 
measurements  using  a  wave  surface  with  2 ad/X =  0.012  and  a  wavelength, 

X  =  2  inches.  These  experiments  provide  the  first  detailed  set  of 
measurements  appropriate  for  testing  models  of  the  wave  induced  variation 
of  the  Reynolds  stress. 

The  analysis  carried  out  by  Thorsness  indicates  that  for  thick 
boundary  layers  the  phase  angle  characterizing  the  shear  stress  variation 
is  a  unique  function  of  a  wave  number,  a  =  2ttv/Xu*,  made  dimensionless 
with  the  kinematic  viscosity,  v,  and  the  friction  velocity,  u  .  The 
amplitude  of  the  shear  stress  variation,  made  dimensionless  using  v  and 
u*,  is  found  to  vary  linearly  with  a  =  au */v.  The  ratio  of  this 
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dimensionless  amplitude  to  a  is  an  unique  function  of  a,  A  major 
limitation  of  the  measurements  of  Thorsness  is  that  they  were  obtained 
in  a  range  of  a  which  was  not  entirely  suitable  to  discriminate  amongst 
various  turbulence  models.  Analysis  suggests  that  the  variation  of  the 
shear  stress  phase  angle  passes  through  a  well  defined  maximum  and 
that  the  prediction  of  this  maximum  should  provide  a  sensitive  test  of 
turbulence  models.  The  experiments  performed  in  this  research  are  an 
improvement  over  the  work  of  Thorsness  in  that  they  covered  a  wide  enough 
range  of  a  to  determine  this  maximum.  This  was  accomplished  by  using 
the  same  wavelength,  solid  waves  of  X  =  2  inches,  and  by  increasing  the 

A 

maximum  value  of  u  by  a  factor  of  four. 

Many  turbulence  models  use  the  kinetic  energy  of  the  turbulence 
fluctuations  as  a  primary  variable  in  estimating  the  Reynolds  stress. 
Very  few  measurements  are  available  of  this  quantity  close  to  a  solid 
surface.  This  lack  of  experimental  data  makes  the  evaluation  of  these 
models  more  difficult.  Therefore,  additional  experiments  were  carried 
out  to  determine  the  streamwise  component  of  the  turbulent  kinetic 
energy  at  the  wave  surface. 

The  measurements  of  the  time  average  shear  stress  and  root  mean 
square  value  of  the  fluctuations  were  obtained  in  a  rectangular  channel 
two  inches  high  and  twenty  four  inches  wide.  The  measurements  were 
obtained  utilizing  an  electrochemical  technique  developed  in  this 
laboratory  by  Reis  [57],  Mitchell  and  Hanratty  [58]  and  Cook  [15]. 

This  technique  makes  use  of  an  electrochemical  solution  as  the  trans¬ 
ducing  medium.  The  solution  flowed  over  a  train  of  eleven  waves 
comprising  the  bottom  wall  of  the  test  section.  Each  wave  had  an 
amplitude  of  0.014  inches  and  a  wavelength  of  2  inches. 
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The  presence  of  waves  on  the  solid  surface  causes  differences  in 
the  turbulence  properties  from  what  would  exist  for  turbulent  flow  over 
a  flat  plate.  This  is  due  to  the  wave  induced  pressure  variations  and 
the  wave  induced  curvature  of  the  streamlines.  A  periodic  pressure 
variation  occurs  along  the  wave  surface  due  to  the  compression  of  the 
streamlines  at  the  wave  crest  and  the  rarefaction  of  the  streamlines 
at  the  wave  trough.  Experimental  studies  of  Jones  and  Launder  [32]  and 
Anderson,  Kays  and  Moffat  [3]  have  shown  that  these  negative  and  positive 
pressure  gradients  lead  respectively  to  a  damping  and  enhancement  of  the 
turbulence.  Similarly  the  alternating  positive  and  negative  curvature 
of  the  streamlines  can  also  lead  to  alternating  enhancement  and  damping 
of  the  turbulence.  The  problem  of  predicting  the  influence  of  a  wavy 
surface  on  the  turbulence  is  further  complicated  in  that  the  turbulence 
does  not  adjust  instantaneously  to  the  change  in  the  pressure  gradient 
and  the  streamline  curvature. 

Two  approaches  are  explored  to  evaluate  the  Reynolds  stress.  The 
first  is  an  extension  of  Thorsness  [70]  Model  D  in  which  he  applied  the 
ideas  of  Loyd,  Moffat  and  Kays  [42]  mixing  length  theory  to  flow  over  waves. 
The  second  is  a  modification  of  the  Jones  and  Launder  [31]  K-e  Model. 

The  advantage  of  the  Loyd  et  al.  mixing  length  approach  is  the  simple 
manner  in  which  the  effects  of  the  pressure  gradient,  streamline  curvature 
and  the  relaxation  can  be  taken  into  account.  The  disadvantage  of  this 
approach  is  the  ad  hoc  manner  in  which  the  relaxation  effects  are 
introduced.  The  K-e  Model  avoids  the  arbitrary  approach  of  introducing 
relaxation  phenomena  by  solving  transport  equations  for  the  turbulence 
properties  which  define  the  turbulent  viscosity. 
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CHAPTER  2 
LITERATURE 

(a)  Effects  of  Pressure  Gradient  and  Curvature 

Experiments  carried  out  in  order  to  understand  the  effects  of 
pressure  gradient  and  streamline  curvature  on  a  turbulent  boundary  layer 
have  led  to  the  development  of  turbulence  models  which  have  been  applied 
to  flow  over  wavy  surfaces.  A  short  review  of  these  experiments,  the 
models  developed,  and  their  applications  is  also  presented. 

An  extensive  review  of  momentum  and  thermal  boundary  layers  subject 
to  pressure  gradients  and  transpiration  has  been  given  by  Kays  and 
Moffat  [33].  Jones  [28],  Launder  and  Stinchcombe  [37],  Badri  Narayanan 
and  Ramjee  [4],  Julien,  Kays  and  Moffat  [32],  Launder  and  Jones  [39], 

Loyd,  Moffat  and  Kays  [42]  and  Jones  and  Launder  [30]  have  measured 
velocity  profiles  and  skin  friction  in  a  turbulent  boundary  layer  which 
has  been  accelerated  by  flowing  through  a  plane-walled  convergent 
channel.  The  interest  in  these  flows  stems  from  the  fact  that  they 
approach  a  state  in  which  the  local  Reynolds  number,  the  skin  friction 
and  the  shape  factor  are  invariant  with  flow  direction  and  therefore 
constitute  one  of  the  simplest  flows  in  which  to  study  the  effect  of 
pressure  gradient.  These  flows  are  a  special  case  of  an  equilibrium 
boundary  layer  and  are  often  referred  to  as  sink  flows  or  asymptotically 

accelerated  boundary  layers.  The  strength  of  the  acceleration  is  charac- 

_2 

terized  by  the  magnitude  of  the  parameter,  K,  defined  as  U  v(dU/dx), 
where  U  is  the  local  free-stream  velocity  and  v  is  the  kinematic  viscosity 
of  the  fluid. 

The  velocity  measurements  indicate  that  for  moderate  acceleration, 

K  -  1  x  10  ^  the  boundary  layer  remains  turbulent;  however,  the  viscous 
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sublayer  becomes  thicker  in  terms  of  the  distance,  y,  made  dimensionless 
with  wall  parameters.  The  velocity  profiles  lie  above  the  universal 
logarithmic  law  of  the  wall.  As  the  acceleration  is  increased  the 
deviation  from  logarithmic  behavior  becomes  more  evident  with  distinction 
between  the  viscous  sublayer  and  the  fully  turbulent  region  of  the  veil,  city 
profiles  becoming  less  clear .  Eventually  a  state  is  reached  in  which  a 
turbulent  boundary  layer  can  no  longer  be  sustained  and  the  flow  is  said 
to  have  undergone  relaminarization. 

Anderson,  Kays  and  Moffat  [3]  have  measured  velocity  profiles  and 
skin  friction  in  unfavorable  pressure  gradients  without  separation.  In 
contrast  to  accelerating  flows,  the  velocity  profiles  remain  logarithmic 
and  experience  a  thinning  of  the  viscous  dominated  sublayer  in  terms  of 
the  dimensionless  distance,  y. 

The  effect  of  streamline  curvature  on  a  turbulent  boundary  layer 

has  been  thoroughly  reviewed  by  Bradshaw  [8]  and  later  by  Gillis  e_t  al. 

[24].  Bradshaw  [7]  has  shown  that  a  flow  subjected  to  mild  longitudinal 

curvature  with  a  ratio  of  boundary  layer  thickness,  6,  to  surface  curvature, 

R  of  6/R  =  1/300  can  significantly  effect  the  length  scale  distribution, 
c’  c 

So  and  Mellor  [66]  have  demonstrated  that  a  convex  surface  inhibits 
turbulence  while  a  concave  surface  enhances  it.  Eskinazi  and  Yeh  [19] 
found  that  the  wall  shear  stress  is  larger  on  the  concave  wall  than  on 
the  convex  wall  of  a  curved  duct.  Ellis  and  Joubert  [18]  have  shown  that 
the  width  of  the  logarithmic  region  is  curvature  dependent.  Convex 
curvature  causes  the  velocity  profiles  to  become  wake  like  at  a  lower 
value  of  y  than  does  concave  curvature.  Ramaprian  and  Shivaprasad  [36,57] 
have  shown  that  the  effects  of  curvature  are  far  more  significant  on 
the  outer  region  of  a  boundary  layer  than  on  the  region  close  to  the  wall. 
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Gillis  et  al.  [24]  carried  out  experiments  to  determine  how  a  boundary 
layer  responds  to  strong  convex  curvature  followed  by  a  flat  section. 

The  experiments  show  that  both  the  shear  stress  in  the  outer  part  of  the 
boundary  layer  and  the  wall  shear  stress  are  strongly  diminished  on 
encountering  the  curved  surface.  When  the  surface  becomes  flat  again 
both  wall  shear  and  shear  stress  profiles  recover  very  slowly  to  flat 
plate  conditions. 

The  first  models  of  turbulent  boundary  layers  close  to  a  solid  wall 
have  used  Van  Driest1 s  [73]  modification  of  Prandtl's  mixing  length 
hypothesis.  This  approach,  which  assumes  a  universality  of  the  wall 
region  in  terms  of  wall  parameters,  fails  to  predict  flows  subject  to 
mild  pressure  gradient  and  streamline  curvature.  Various  workers  have 
proposed  modifications  to  the  Van  Driest  formula  to  better  account  for 
these  effects. 

Patankar  and  Spalding  [53]  proposed  that  the  local  value  of  the 
shear  stress  rather  than  the  wall  value  be  used  in  the  exponent  of  the 
Van  Driest  damping  function.  This  modification  has  the  correct  quali¬ 
tative  behavior.  In  favorable  pressure  gradients  the  shear  stress 
decreases  from  its  value  at  the  wall,  so  that  this  formulation  does 
result  in  a  thickening  of  the  viscous  wall  region.  However,  calculations 
carried  out  using  this  formulation  show  that  it  does  not  produce  a 
large  enough  effect. 

Launder  and  Jones  [38]  ,  Cebeci  and  Smith  [12],  Julien  et  al.  [32]  and 
Loyd  etal.  [42]  began  experimenting  with  the  idea  that  the  constant  A  in  the 
Van  Driest  mixing  length  model  is  related  to  the  thickness  of  the  viscous 
sublayer  in  wall  coordinates  which  in  turn  depends  on  the  dimensionless 
pressure  gradient,-^-  .  These  workers  deduced  a  functional  dependency 
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of  A  on  ^  by  examining  a  large  number  of  equilibrium  velocity  profiles 
dx 

subject  to  both  positive  and  negative  pressure  gradients. 

In  nondimensional  flows  where  the  pressure  gradient  is  changing 
rapidly  Loyd  et  al.  [42],  Julien  et  al .  [32]  and  Launder  and  Jones  [38] 
have  proposed  that  a  first  order  lag  equation  be  used  to  predict  an 
effective  sublayer  thickness,  Aeff. 

Bradshaw  [7]  has  modified  the  Van  Driest  mixing  length  proposal  to 
account  for  the  effect  of  streamline  curvature  on  the  turbulence.  The 
flat  wall  mixing  length  is  multiplied  by  an  empirical  function  built 
around  the  curvature  Richardson  number.  In  a  situation  in  which  the 
surface  curvature  is  not  constant,  Bradshaw  suggests  that  an  effective  radius 
of  curvature  be  computed  from  a  first  order  lag  equation. 

Jones  and  Launder  [29,31]  abandoned  the  mixing  length  approach  by 
arguing  that  the  prediction  of  nondimensional  flows  could  not  be  achieved 
with  a  transport  hypothesis  based  so  firmly  on  equilibrium  notions.  They 
assumed  that  the  turbulent  viscosity  is  the  product  of  the  square  root  of 
the  turbulent  kinetic  energy  and  an  appropriate  length  scale.  The  turbu¬ 
lence  length  scale  is  calculated  by  solving  transport  equations  for  the 
turbulent  kinetic  energy  and  the  turbulent  dissipation  rate.  Launder 
et  al .  [40]  have  modified  this  model  to  account  for  streamline  curvature, 
by  redefining  the  curvature  Richardson  number  in  terms  of  turbulence 
quantities.  Several  workers  have  developed  alternate  higher  order  closure 
schemes.  These  have  been  reviewed  by  Reynolds  [59,60] 

(b)  Models  for  Flow  Over  Wavy  Surfaces 

Benjamin  [6]  and  Miles  [46]  considered  the  wave  induced  flow  caused 
by  small  amplitude  waves.  They  formulated  the  problem  in  a  curvilinear 
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coordinate  system  and  included  the  effects  of  turbulence  only  in  the 
specification  of  the  mean  velocity  profile. 

Hussain  and  Reynolds  [27]  incorporated  the  effects  of  turbulence 
on  the  perturbed  flow.  The  eddy  viscosity  distribution  is  assumed  to 
be  set  up  by  the  mean  flow.  The  mean  eddy  viscosity  is  then  assumed 
to  act  on  the  perturbed  velocity  gradient  generating  a  wave  induced 
Reynolds  stress. 

Other  workers  to  examine  flow  over  wavy  surfaces  include  Davis  [16] 
and  Townsend  [72].  Their  models  have  been  reviewed  by  Thorsness  [70], 

Thorsness  [70]  formulated  the  problem  in  a  curvilinear  coordinate 
system  and  investigated  several  turbulence  models  for  the  wave  induced 
Reynolds  stress.  Three  particular  models  are  of  interest.  They  are 
labeled  as  Model  A,  Model  C,  and  Model  D.  Model  A  is  essentially  the 
quasilaminar  model  of  Benjamin  described  above.  Model  C  evaluates  the 
wave  induced  Reynolds  stress  by  using  the  Van  Driest  mixing  length  model. 
A  wave  induced  eddy  viscosity  results  because  the  shear  stress  used  in 
the  damping  function  is  the  local  wall  shear  stress.  Model  D  is  an 
adaptation  of  the  mixing  length  model  of  Loyd  et  al.  This  model  uses 
the  wave  induced  pressure  gradient  and  wall  shear  stress  in  evaluating 
the  mixing  length. 

Markatos  [44]  solved  the  full  nonlinear  problem  using  a  curvilinear 
coordinate  system.  In  order  to  avoid  the  difficulty  of  modeling  the 
turbulence  in  the  viscous  wall  region  the  high  Reynolds  number  form  of 
'the  Jones  and  Launder  two  equation  model  is  matched  to  the  logarithmic 
law  of  the  wall. 

Cary  elt  al.  [11]  have  used  the  mixing  length  model  of  Loyd  et  al. 
and  the  two  equation  model  of  Jones  and  Launder.  Their  study  is 
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primarily  concerned  with  the  investigation  of  a  wavy  surface  as  a 
possible  drag  reducing  device.  Their  analysis  suggests  that  neither 
of  the  two  models  describe  flow  over  large  amplitude  waves. 

Caponi  er  al.  [10]  have  developed  a  model  for  laminar  flow  over 
an  arbitrarily  shaped  periodic  surface.  An  orthogonal  transformation 
is  used  to  map  the  physical  domain  under  consideration  to  a  rectangular 
region.  Because  of  the  periodic  nature  of  the  problem,  the  dependent 
variables  are  expanded  in  terms  of  Fourier  series.  The  model  sucess- 
fully  predicts  viscous  flow  over  both  a  moving  liquid  and  solid  wave 
surface. 

Mclean  [45]  has  extended  the  laminar  flow  calculations  of  Caponi 
et  al.  to  include  turbulent  flow  by  using  the  mixing  length  model  of 
Loyd  _et  al.  A  comparison  of  this  theory  with  the  shear  stress  measure¬ 
ments  of  Thorsness  and  Zilker  indicates  that  the  model  successfully 
predicts  flow  over  small  amplitude  waves.  However, in  the  case  of 
large  amplitude  waves  the  discrepancy  between  the  theory  and  experiments 
is  more  apparent. 
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CHAPTER  3 
THEORY 

In  this  chapter  a  theoretical  framework  for  the  description  of 
turbulent  flow  over  a  small  amplitude  wave  surface  is  presented.  The 
problem  is  formulated  in  a  boundary  layer  coordinate  system.  A  number 
of  models  for  the  wave  induced  Reynolds  stresses  are  developed.  The 
results  of  the  alternate  approaches  and  their  comparative  success  in 
predicting  the  experimental  results  is  reserved  for  presentation  in 
Chapter  6  and  Chapter  7.  The  various  numerical  techniques  used  in 
integrating  the  equations  are  presented  in  Chapter  5. 

(a)  Coordinate  System  and  Basic  Equations 

The  boundary  layer  coordinate  system  used  in  formulating  the 
problem  is  shown  in  Figure  3.1.  The  x  direction  is  taken  parallel  to 
the  wave  surface  while  the  y  axis  is  perpendicular  to  it.  For  con¬ 
venience  x  =  0  is  taken  as  the  wave  crest  with  the  positive  x  axis 
being  in  the  direction  of  flow.  The  flow  field  is  assumed  to  be  two 
dimensional. 

Unless  otherwise  stated  all  the  variables  are  made  nondimensional 

with  respect  to  wall  parameters.  Velocities  are  made  dimensionless 

*  rzi - - — 

with  respect  to  the  friction  velocity,  u  =  /t  / p  .  Lengths  are 

wd 

made  dimensionless  with  respect  to  v/u  ,  where  v  is  the  kinematic 

viscosity.  Pressure,  stresses  and  the  turbulent  kinetic  energy  are 

*2 

considered  multiples  of  p  u 

The  time  averaged  continuity  and  momentum  equations  in  the  boundary 
layer  coordinate  system  with  surface  curvature,  k,  are: 
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where  U  and  V  are  the  velocities  in  the  x  and  y  directions  respectively. 


The  quantities,  -u'u'  ,  -u'v'  and  -v'v'  are  the  Reynolds  stresses  and  P 


is  the  static  pressure. 
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(b)  Linearized  Equations 

The  necessary  conditions  for  the  flow  field  to  be  described  by  a 
linear  form  of  the  Navier  Stokes  equations  are  that  the  wave  amplitude 
be  small  compared  to  both  the  wavelength  and  the  boundary  layer 
thickness.  These  conditions  are  developed  in  Appendix  A. 

The  time  averaged  velocities,  U,  V,  the  turbulent  stresses,  ,  and 
the  pressure  gradients,  3P/3x,  3P/3y,  are  assumed  to  be  the  sum  of  a 
component  averaged  over  a  wavelength  at  a  constant  value  of  y  and  a 
periodic  spatially  varying  wave  induced  component.  The  general  form 
of  the  wave  induced  component  is  anelax  where  n  is  a  complex  number 
whose  real  and  imaginary  parts  are  at  most  functions  of  y.  Thus  the 
velocities,  stresses  and  pressure  gradients  are  given  by 


U  =  U(y)  +  au(y)e1QlX 


(3.4.1) 


„  .  iax 

V  =  av(y)e 


(3.4.2) 


(3.4.3) 


■dP.(y)  eiotx 
dy 


ap 

3y 


(3.4.4) 


(3.4.5) 


If  the  equations  (3.4.1)-(3.4.5)  are  substituted  into  the  continuity  and 

2 

momentum  equations  (3.1) -(3. 3)  and  the  terms  of  0(a  )  are  neglected,  a 
system  of  equations  which  is  linear  in  the  wave  induced  components  results. 
The  continuity,  the  x  momentum  and  the  y  momentum  equations  are  respectively 


r  .  -  ,  n  ictx  n 
a  I  iau  +  v  J  e  =*  0  , 


(3.5) 
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r  »  tta  i  att  i  i  10X 
a  [  iaUu  +  vU  ]  e  - 


+  U"  +  R' 
8x  xy 


+  a  [  -iap  -  a^u  +  Uf  +  a^y  +  u,f 


8x 


+  2R  +  iar  +  rT  ]  e^aX 

xy  xx  xy 


(3.6) 


r  .  77-  2  —2  iax 

a  [  iaUv  -  a  U  ]  e 


-f  +  R’ 

3y  yy 


+  a  [  v"  -  a^v  -  ia"^  U-p'  +  r*  +  R  a^  +  iar 

yy  yy  xy 


2  -5-  ,  iax 

a  R  ]  e 

xx 


(3.7) 


The  primes  denote  differentiation  with  respect  to  the  y  direction.  The 
continuity  equation  (3.5)  can  be  satisfied  by  the  introduction  of  the 
stream  function 


¥  * 


U  (y)  dy  +  a  F(y)  e10lX  , 
'o 


with 


8^ 

8y 


(3.8) 


(3.9.1) 


and 


V 


zi  II 

h  8x 
x 


(3.9.2) 
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where  h  and  h  are  the  metric  functions.  The  linearized  metrics 
x  y 

are: 


h  =  1  (3.10.1) 

x 

and 

h  =  1  +  aa2yelaX-  (3.10.2) 

y 

Substituting  equations  (3.10.1)  and  (3.10.2)  into  equations  (3.9.1)- 
(3.9.2)  and  using  the  definition  of  the  stream  function  equation  (3.8), 
the  U  and  V  velocities  to  0  (a)  are  respectively, 

U=U+aF'eiaX  (3.11.1) 

and 

V  =  -aictFeiotX  .  (3.11.2) 

The  linearized  momentum  equations  (3.6)  and  (3.7)  can  be  combined 

to  eliminate  the  pressure  terms  in  the  following  manner.  Equation  (3.6) 

2  iotx 

is  multiplied  by  (1  +  aa  y  e  )  and  then  differentiated  with  respect 
to  y.  Equation  (3.7)  is  differentiated  with  respect  to  x  and  multiplied 
by  (-1).  The  resulting  equations  are  then  added.  Substituting  the 
definitions  of  the  velocities  in  terms  of  the  stream  function  (3.11.1) 
and  (3.11.2)  into  the  resulting  equation  and  collecting  terms  of  similar 
order  the  following  equations  defining  F(y)  are  obtained: 

U™  +  R"  =  0  ,  ‘  (3.12) 

xy 
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ict  [  U (F"  -  ct2F)  -  U"F  +  a2  U2  ] 


=  F""  -  2a2F"  +  a4F  +  2a2  U"  -  a4  U  +  R 


(3.13) 


R =  ia3  R 


3  - 
-  ia  R 


(3.14) 


XX 


yy 


Equation  (3.13)  is  similar  to  the  well  known  Or r  —  Sommer f eld  equation 
used  in  stability  calculations  except  for  the  appearance  of  additional 
terms  which  arise  due  to  the  use  of  the  curvilinear  coordinate  system 
and  due  to  the  inclusion  of  the  Reynolds  stresses. 

Solution  of  equations  (3.12)  and  (3.13)  requires  the  specification 
of  the  Reynolds  stresses  and  the  boundary  conditions.  The  average 
velocity  profile,  U(y),  and  the  average  Reynolds  stress,  R^(y),  appearing 
in  equations  (3.12)  and  (3.13)  are  taken  to  be  the  same  as  would  exist  if  the 
surface  were  flat.  This  specification  is  consistent  with  the  linearization 
assumption  in  that  the  leading  order  wave  averaged  terms  are  the  same  in 
both  a  boundary  layer  and  cartesian  coordinate  system.  Furthermore  it 
is  assumed  that  the  wave  averaged  shear  stress,  x(y),  does  not  vary  in  the 
y  direction.  This  implies  that  the  wave  averaged  pressure  gradient, 

3P/3x  is  zero.  Equation  (3.12)  can  therefore  be  integrated  once  to  yield 


U"  +  R* 


(3.15) 


0 


xy 


The  specification  of  the  Reynolds  stress  is  dealt  with  in  the  next  section. 

At  the  wave  surface  no  slip  and  no  penetration  boundary  conditions 
are  invoked. 


F  =  0  and  F!  =  0  at  y  =  0 


<  3 .  lb ,) 
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The  outer  boundary  condition  is  such  that  the  flow  field  far  away  from 
the  surface  is  unaffected  by  the  presence  of  the  wave.  Therefore  the 
velocity  field  in  the  boundary  layer  coordinate  system  must  be  equated 
to  an  undisturbed  velocity  field  that  would  exist  if  the  wave  surface 
were  not  present.  This  is  accomplished  in  two  steps.  Firstly  the 
velocity  field  in  the  boundary  layer  coordinate  system  is  transformed 
into  cartesian  coordinates.  Secondly  the  transformed  velocity  field 
is  then  equated  to  a  linearized  velocity  field  in  cartesian  coordinates, 
see  Appendix  B.  The  linearized  boundary  conditions  for  the  flow  far 
from  the  wave  surface  are  then  obtained  as 

F  =  U  and  F,==  U*  for  large  y.  (3.17) 


The  shear  stress  at  the  surface  is  evaluated  from  the  component  of 

the  rate  of  strain  tensor,  S  , 

xy 


T 

w 


2  S 


xy 


h  A 

fu] 

1  3V " 

_hx  3y 

W 

h  3x 
x  J 

(3.18) 


2 

Neglecting  terms  of  0(a  )  and  using  the  surface  boundary  conditions 
equation  (3.16), 


x  =  x (0)  +  ax (0)elaX  =  U'(0)  +  aF"(0)eiax  .  (3.19) 

W 


The  wave  induced  variation  of  the  wall  shear  stress  is 

ax(0)eiaX  =  aFn(0)e^aX  .  (3.20) 

The  wave  induced  variation  of  the  pressure  at  the  solid  surface  is 
found  by  evaluating  equation  (3.6)  at  y  =  0 
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ap(0)elax  =  ^  [F,M  (0)  +  a2U'  (0)]  aelaX 


(3.21) 


(c)  Models  for  Turbulent  Stresses 

c. 1  Constant  Reynolds  Stress  Assumption  Model  A 

The  simplest  assumption  to  make  is  the  quasi-laminar  assumption  of 
Miles  and  Benjamin  that  the  effects  of  turbulence  enter  the  problem  only 
in  the  specification  of  the  mean  velocity  profile  U(y).  In  this  case 
ft  =0.  A  slightly  different  approach  is  taken  here  in  that  only  the 
fluctuating  component  of  the  Reynolds  stress,  r^  ,  is  set  equal  to  zero. 
Consequently 


R  =  3a2  R  ' 


xy  9 


(3.22) 


if  the  normal  stresses,  R  and  R  ,  are  neglected. 

’  xx  yy 

For  a  thick  boundary  layer 


UM  »  -Rf 
xy 


(3.23) 


so 


R  =  -3a2  U  " 


(3.24) 


Thorsness  [70]  has  argued  that  the  constant  stress  model  formulated 
in  a  boundary  layer  coordinate  system,  which  he  designated  as  Model  A, 
is  closely  related  to  the  assumption  that  the  Reynolds  stresses  are 
frozen  along  a  streamline  as  discussed  by  Davis  [16].  In  a  cartesian 
coordinate  system  the  constant  stress  assumption  reduces  to  the  quasi- 
laminar  model  of  Miles  and  Benjamin. 
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A  major  shortcoming  of  the  constant  stress  model  (Model  A)  is 
that  there  is  no  interaction  between  the  wave  and  the  turbulence.  Several 
models  are  now  explored  which  take  this  interaction  into  account.  These 
models  are  all  based  on  the  assumption  that  the  Reynolds  stresses  are 
described  by  a  Newtonian  constitutive  equation  of  the  form 


1  2  . 
t  q 


(3.25) 


Z 

where  q  is  the  turbulent  kinetic  energy  and  v  is  a  turbulent  or  eddy 

viscosity  and  S^.  is  the  rate  of  strain  tensor  defined  in  Appendix  C. 

Since  turbulent  stresses  are  important  only  close  to  the  wave 

surface,  use  is  made  of  the  boundary  layer  assumption  that  the  normal 

stresses, R  and  R  , can  be  neglected.  Here  the  only  component  con- 
xx  yy 

tributing  to  the  Reynolds  stress  tensor  is  R  where 

xy 


R  =  —  2S  .  (3.26) 

xy  v  xy 


In  order  to  account  for  the  spatial  variation  of  the  turbulence  properties 
along  the  wave  surface  the  concept  of  a  wave  induced  eddy  viscosity  is 
introduced 


v 


t  lax 
a  —  e 
v 


(3.27) 


Decomposing  the  rate  of  strain  S 
components 


into  its  wave  averaged  and  fluctuating 


iax 

e 


S  *  S  +  a  s 
xy  xy  xy 


(3.28) 


21 


the  Reynolds  stress,  R  ,  is  related  to  the  eddy  viscosity  and  the  rate 
of  strain  as  follows 


R  =  —  2  S  +  a 

xy  v  xy 


vt  -  Vt  - 

—  2  S  +  —  2  s 

v  xy  v  xy 


iax 


(3.29) 


The  wave  induced  variation  of  the  Reynolds  stress  is 


-  iax 
are  =  a 
xy 


(  v 


—  2  S  +  —  2  s 
v  xy 


\ 

V 


xy 


iax 


(3.30) 


The  average  and  wave  induced  rates  of  strain  are  respectively 


2  S  =  U' 
xy 


(3.31) 


o  -  iax  2  2  — .  iax 

2a  s  e  =  a  (F  +  a  F-aU)e 
xy 


(3.32) 


The  Reynolds  stress  term,/?,  equation  (3-14),  can  now  be  written  as 


/?=  3a 


r  —  > 

Vt  771 

r 

,  2 

—  U1 

+  a 

V 

s  4 

Vt  2  2-  Vt  - 

—  (F"  +  a  F  -  a  U  )  +  —  U' 
V  V 


Vt  2  2-  Vt  - 

—  (F"  +  a  F  -  a  U)  +  —  U* 

V  V 


(3.33) 


Before  equation  (3.12)  to  (3.14)  can  be  solved  for  the  fluctuating  flow 
the  eddy  viscosity  must  be  specified.  Two  approaches  are  explored. 

c. 2  Zero  Equation  Models 

At  this  level  of  turbulence  modeling  the  Prandtl  mixing  length 
hypothesis  is  frequently  used. 
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v 


t 


2  S 


xy 


(3.34) 


where  l  is  "the  mixing  length."  Close  to  a  solid  boundary  where 
o 

viscous  effects  are  important  several  workers  have  adopted  the  Van  Driest 
mixing  length  model 

Aq  *  <y  [1  “  exP  (_Dm)]  (3.35) 

where  k  is  the  Von  Karman  constant  and,  D^jisa  damping  function  introduced 
to  account  for  the  effects  of  the  wall  in  suppressing  turbulent  transport. 

Van  Driest 's  original  proposal  was  that  Dm  =  ^  Tw  where  a  value  of 
A  =  A  =  26  is  suggested  for  flow  over  a  flat  plate.  Thorsness  [70]  has 
applied  this  formulation  of  the  Van  Driest  mixing  length  hypothesis,  which 
he  labeled  as  Model  C,  to  flow  over  small  amplitude  waves.  This  formu- 
lation  is  not  expected  to  be  correct  for  flow  over  wavy  surfaces  because 
the  presence  of  pressure  gradients  causes  a  large  variation  of  x  with 
distance  from  the  boundary  and  a  drastic  change  in  the  production  of 
turbulence  in  the  wall  region.  This  can  be  taken  into  account  by  redefin¬ 
ing  the  damping  function  as 


-'m 


=  ^T1/2 


(y) 


(3.36) 


where  the  local  shear  stress,  x(y),  is  used  instead  of  the  value  at  the 

surface,!  .  The  local  shear  stress  is  the  sum  of  viscous  and  turbulent 
w 

stresses 


S 

xy 


T 


(3.37) 
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For  small  amplitude  waves  the  local  shear  stress  can  be  separated 
into  an  average  and  a  wave  induced  component.  The  average  shear  stress 
is  simply  equal  to  its  value  at  the  surface 


T  = 


1  +  — 

V 


2  S  =  1 
xy 


and  the  wave  induced  component  of  x  is 


(3.38) 


*  *  \  iax  - 

at(y)  e  =  2a 

+ 

vti 

i  +  — 

s 

v  xy 

v 

xy 

The  mixing  length  to  0(a)  is 


(3.39) 


i 

o 


i 

o 


+  a  Z 

o 


iax 

e 


(3.40) 


Using  the  definitions  of  vt>  Iq  and  Dm,  equations  (3.34),  (3.35)  and 
(3.36)  respectively,  and  neglecting  terms  of  O(a^)  the  average  components 
of  the  mixing  length,  and  the  eddy  viscosity,  v^,  are  respectively 

Iq  =  Ky  [1  -  exp  (-y/A)]  (3.41) 


and 


v.  =  l  2  U' 
t  o 


(3.42) 


The  wave  induced  components  of  Zq  and  v  are  respectively 


*  ictx  _  y 

a  l  e  =  a<y  exp  (-y/A)  — 


12  '  ^ 


(3.43) 


and 
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aO  . 
t  iax 

_  A 

s  — 

xy  2  exp  (— y/  A) 

r  y 

X 

X 

vt  iax 

_  D 

- e  =  a 

V 

Sxy  [1  -  exp  (-y/A)] 

.  A  1 

2 

V 

(3.44) 


where 


JDL  = 


F"  + 


2 

a  F  - 


a2  U 


S 

xy 


U' 


(3.45) 


Eliminating  the  eddy  viscosity,  between  equations  (3.39)  and  (3.44) 
an  explicit  result  for  the  variation  of  the  wave  induced  shear  stress 
is  obtained. 


at (y)e 


V  v  ' 

a 

2—  +  1 

2  s 

iax 

A  vt  J 

xy  J 

—  2  S  exp  (-y/A)  y 
v xy  r 

[1  -  exp  (-y/A) ]  A 


iax 

e 


(3.46) 


Equation  (3.44)  and  (3.46)  defining  v  and  x  respectively  are  designated 
as  Model  C  .  If  t  is  evaluated  at  the  wall  Model  C  reduces  to  Model  C. 

For  equilibrium  boundary  layers  Loyd  et:  al.  suggest  the  following 
functional  dependence  of  A  on  pressure  gradient 


A  -  A 


1  +  k. 


dP 


dx 


+  k„ 


dP 


dx 
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(3.47) 


where  A  is  the  flat  plate  value. 

For  flow  over  small  amplitude  waves  A  can  be  expressed  as 


A  =  A  +  aAe 


iax 


(3.48.1) 
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where 


aAe^aX  =  aA  iap(o)e‘'"aX  . 
The  damping  function  is 


(3.48.2) 


r,  ,  -  iax  1/2 

d.  -  y  [1/-  aTe  1 


A  +  a  A  e 


iax 


(3.49) 


2 

Using  equation  (3.35)  and  (3.36)  and  neglecting  terms  of  0(a  )  the  wave 
induced  mixing  length  and  eddy  viscosity  are  respectively 


ai  e 
o 


^ax  =  a<y  exp  (-y/A) 

A 


_A 

A 


lax 


and 


(3. 50) 


°t  iax 
a  —  e 
v 


s 

S 

xy 


2  exp  (-y/A)  fa 
[1  -  exp  (-y/5)]  (A 


iax 

e 


(3.51) 


Eliminating  the  eddy  viscosity  between  equations  (3.39)  and  (3.51)  the 
wave  induced  shear  stress  is 


ax (y)e 


iax 


2  —  +  1 

v 


2  s 


xy 


Vt  - 

2  —  U'  exp  (-y/A) 
[1  -  exp  (-y/A)  ] 


iax 

-  e 


(3.52) 


r  !t 

v 


1  - 


UT  exp  (-y/A) 

f  ' 

y 

I 

-  exp  (-y/A) ] 

A 

V.  > 

Equations  (3.51),  (3.52  and  (3.48.2)  represent  an  equilibrium  eddy 
viscosity  model  for  flow  over  a  small  amplitude  wave. 
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In  situations  where  the  pressure  is  varying  rapidly  Loyd  et_  al. 
[42]  (see  literature  survey)  have  suggested  that  the  flow  close  to  the 
wall  sees  an  effective  pressure  gradient, 


_d_ 

dx 


dP 


dx 


eff 


dP 

dx 


dP 

dxj 


eff 


k 


L 


(3.53) 


where  k  is  a  lag  constant. 

JL< 

Solving  for  a  periodic  pressure  gradient 


-  iax 

~  iax  ape 
ap  e  =  - - 

ett  [1  +  iakj 


and  thus  replacing  p  with  Pe££  equation  (3.48.2)  becomes 


(3.54) 


a  A  r  -  e 
eff 


iax 


a  A  k^i  a  p 

[1  +  i  a  k  ] 
L 


(3.55) 


Thorsness  [70]  has  labeled  the  above  formulation  for  A^^  together  with 
equation  (3.51)  and  equation  (3.52)  evaluated  at  the  wall  Model  D.  Abrams, 
Frederick  and  Hanratty  [2],  using  the  local  value  of  the  shear  stress 
equation  (3.52)  together  with  equation  (3.51)  and  (3.55),  have  labeled  this 
formulation  Model  D  . 

A  slightly  different  approach  is  also  explored  in  this  research. 
Instead  of  accounting  for  the  nonequilibrium  effects  in  the  boundary 
layer  by  relaxing  the  pressure  gradient,  the  complete  Reynolds  stress  term, 
R ,  is  relaxed,  by  introducing  an  effective,  ^e-£»  satisfying  the  following 
first  order  rate  equation. 
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d  R  iax 

—  n  a  e 
dx  ef  f 


<* 


iax 
a  e 


(3.56) 


Solving  for  a  periodic  variation  yields, 


R 

eff 


R 

[1  +  i a  kR] 


(3.57) 


Equations  (3.52),  (3.51),  (3.48.2)  and  (3.57)  for  vt»  t  ,  A,  and  R  define 
a  nonequilibrium  relaxation  model  for  flow  over  a  small  amplitude  wave. 

The  mean  velocity  profile  to  be  used  with  the  above  turbulence 
models  is  obtained  by  integrating  the  wave  averaged  momentum  equation 
(3.15).  This  results  in  the  well  known  Van  Driest  velocity  profile. 


U(y) 


'7  2  dy 

4>  1  +  /  (1  +  4<2y^[l  -  exp  (-y/A)]2) 


(3.58) 


c.3  K-e  Model 

The  basic  high  Reynolds  number  K-e  model  used  in  this  thesis  is 
the  model  proposed  by  Jones  and  Launder  [29,31].  Several  workers  have 
extended  the  basic  model  to  include  the  low  Reynolds  number  viscous  wall 
region.  Patel  £t  al.  [54]  have  tested  many  of  these  models  against  data 
sets  obtained  in  a  variety  of  external  pressure  gradients.  The  most 
successful  low  Reynolds  number  model  tested  is  the  one  proposed  by 
Chien  [13].  This  is  the  model  adopted  for  use  in  this  thesis. 

The  K*-e  model  requires  that  a  partial  differential  equation  be 
developed  for  the  turbulent  kinetic  energy 


(3.59) 
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and  for  the  isotropic  turbulent  energy  dissipation  rate 


e  = 


3ul  3u! 
_ 1  _ 1 

3x .  3x. 
J  3 


(3.60) 


Equations  for  the  kinetic  energy  and  dissipation  rate  are  presented  in 
Appendix  C  together  with  the  necessary  closure  assumptions.  The  closed 
form  of  the  energy  and  dissipation  equations  in  boundary  layer  coordinates 
are  respectively 
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where 


f  =  1-0. 4/1. 8  exp  (-K^/36e^) 
f d  =  exp  ^-C4 y Tw> 


Terms  I 
Terms  II 
Terms  III 

Terms  IV 
Terms  V 


advection  of  K  or  e  by  time  averaged  flow 
viscous  diffusion  of  K  or  e 
turbulent  diffusion  of  K  or  e  by  pressure 
and  velocity  fluctuations 

production  of  K  or  e  by  time  averaged  flow 
dissipation  of  K  or  e 


(3.63-1) 

(3.63.2) 


The  eddy  viscosity  v  is  defined  as 


2 

v  C 

—  -  — -  [1  -  exp  (— C~  y  x  )  ]  . 

v  e  3  w 


(3.64) 


In  summary,  the  turbulence  model  is  assumed  to  be  governed  by  equations 
(3.61)-  (3.64)  with  the  constants  a ^  *  1,  0£  =  1.3,  =  0.09,  *  1.35, 

*  1.8,  =  0.0115,  and  =  0.5.  The  values  used  are  those  suggested 

by  Chien  [13]. 
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For  flow  over  small  amplitude  waves  it  is  possible  to  define  a 
mean  and  fluctuating  component  of  the  turbulent  energy  and  dissipation 
rate 


K  =  K(y)  +  ak(y)e^aX  (3.65) 

e  =  e(y)  +  ae(y)eiax  .  (3.66) 

Substituting  equations  (3.65)  and  (3.66)  into  equations  (3 . 60)- (3 . 61) , 

2 

collecting  terms  of  like  order  and  neglecting  0(a  )  terms,  the  following 
equations  for  the  mean  and  fluctuating  energy  and  dissipation  rate 
result. 

Mean  energy  equation: 

(3.67) 

(3.68) 

where 


[(1  +  vt/voR)  K'*j 


/vO  K'  I  +  R  U *-  e  - 

K7  |  xy  y2 


Mean  dissipation  equation: 
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. 36e2  - 

fd  =  exp  (-C^y) 


\  K2 

—  =  C  [  1  -  exp  (-C,y)  ] 

l1  e  J 


(3.69.1) 


(3.69.2) 


v 


(3.69.3) 
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are 


K=e  =  0  at  y  =  0 


(3.70) 


and  at  large  y  the  energy  and  dissipation  distribution  must  be  such  as 
to  yield  a  constant  stress  region 

-1/2 


K  =  C 

y 

e  =  l/<y  . 

Fluctuating  energy  equation : 


for  large  y.  (3.71) 
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Fluctuating  dissipation  equation  : 
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where 


*  —  2 
vfc  C  Kz 


4-  [  1  -  exp  (-C_y)  ] 
e  J 


r  2  k  _  £  ' 
.  K  l  4 


K2 


+  c3c  exP  (~c3y)  f  F"(0) 


(3.74) 


and 


P,  °  2  R  s  +  2  r  S 
k  xy  xy  xy  xy 


+  R  s  +  R  s 
xx  xx  yy  yy 


(3.75) 


Assuming  that  the  normal  stresses  are  small  and  using  equation 
(3.30)  the  production  term  becomes 


P 


k 


(3.76) 


The  velocity  field  is  obtained  by  solving  the  stream  function  equation 
(3.13)  simultaneously  with  equations  (3. 71)-(3. 76) . 
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The  boundary  conditions  at  the  wave  surface  and  at  large  values  of  y 
for  the  fluctuating  energy  and  dissipation  are 

k  =  e  =  0  at  y  =  0  (3.77.1) 

and 


k  =  e  =  0  for  large  y  (3.77.2) 

Since  equations  (3.72)  and  (3.73)  are  singular  at  y  =  0  the  boundary 
conditions  at  the  wave  surface  cannot  be  used  in  the  form  presented  in 
equation  (3.77.2),  If  a  Taylor  series  expansion  is  made  on  equations  (3.70) 
and  (3.71)  about  y  «  0  the  following  boundary  conditions  are  arrived  at 


at  y  =  6d 


(3.78.1) 


0 


(3.78.2) 


where  6 


d 


is  a  small  distance  from  the  wave  surface. 


(d)  Streamline  Curvature 

The  influence  of  streamline  curvature  is  taken  into  account  at  the 
zero  equation  level  of  turbulence  modeling  by  modifying  the  plane  shear 
mixing  length  relation  equation  (3.35)  in  the  manner  suggested  by  Bradshaw. 

I  =  %  (1  -  3  R.  )  (3.79) 

o  Cl 

c 


Here  g^  is  an  empirical  constant  and  IL  is  the  curvature  Richardson 
number  defined  as 


R. 
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2U 
3U 
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(3.80) 
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where  1/R  is  the  curvature  of  the  streamlines.  For  small  amplitude  waves 


1/R 


c 


a 


(F  -  U  )  a2eiax 
U 


(3.81) 


in  boundary  layer  coordinates.  The  details  of  this  derivation  are 
given  in  Appendix  D. 

Bradshaw  suggests  that  an  effective  radius  of  curvature, 
satisfying  a  first  order  lag  equation  be  used 


_d _ 1 

dx  R 


c 


eff 


R 

c 


1 


R 

cef  f 


k 

c 


(3.82) 


Here  is  the  curvature  lag  constant.  Solving  the  above  equation  for 
a  periodic  variation,  the  effective  radius  of  curvature  of  the  streamlines 
is 


1  _  1 
R  R  (1  +  iak  ) 

cef  f  c 


(3.83) 


Using,  1/R  in  equation  (3.80)  and  neglecting  terms  of  0(a  ),  the 
cef  f 

curvature  Richardson  number  can  be  written  as 
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(3.85) 


35 

The  wave  induced  mixing  length,  eddy  viscosity  and  shear  stress 
distribution  incorporating  the  effect  of  curvature  are  respectively 
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(e)  Finite  Boundary  Layer  Calculation 

The  calculations  thus  far  have  assumed  that  the  boundary  layer  is 
sufficiently  deep  such  that  in  the  region  where  the  wave  induced  flow 
is  negligible  the  average  velocity  profile  is  still  logarithmic.  The 
problem  has  therefore  been  characterized  by  a  single  parameter  the  wave 
number, a.  In  this  section  the  requirement  of  an  infinite  boundary  layer 
thickness  is  relaxed. 

e. 1  Mean  Flow 

The  average  velocity  profile  is  assumed  to  be  given  by  the  Coles 
proposal  (see  Cebeci  and  Smith  [12]). 
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2n 


sin 


(3.89) 


where  II  =  0.55  for  a  zero  pressure  gradient  boundary  layer  and  6  is  the 
boundary  layer  thickness. 

e. 2  Eddy  Viscosity 

A  finite  boundary  layer  consists  of  two  parts:  a  fully  turbulent 
inner  region,  where  the  Reynolds  stress  distribution  is  calculated  from 
equations  (3.51)  and  (3.52)  for  the  average  and  wave  induced  flows 
respectively  and  an  outer  region,  where  Cebeci  [12]  suggests  that  the 
mixing  length  is  given  as 

Z  =  y  6  y  £  y  £  6  •  (3.90) 


Here  y  is  an  empirically  determined  constant  approximately  equal  to 
0.075  and  y is  obtained  from  the  continuity  of  £o.  The  eddy  viscosity 
is  found  using  (3.34)  to  be 
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Neglecting  terms  of  0(a  )  the  average  eddy  viscosity  is 


=  y262  U' 
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(3.92) 


and  the  wave  induced  eddy  viscosity  is 
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Using  equation  (3.30)  the  wave  induced  Reynolds  stress  in  the  outer 
region  of  the  boundary  layer  becomes 


,  a  v 

*  iax  t  ,  A  lax 

are  =  -  4  s  e 

xy  v  xy 


y  <  y  <  5 
;c 


(3, 


94) 


The  details  of  the  finite  boundary  layer  calculation  are  presented  in 
Appendix  E. 
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CHAPTER  4 

EXPERIMENTAL  MEASUREMENTS  AND  RESULTS 

The  purpose  of  the  experimental  effort  was  to  obtain  data  which 
can  be  used  to  compare  the  various  turbulence  models  developed  in 
Chapter  3  for  the  wave  induced  Reynolds  stresses.  The  data  consist  of 
time  average  shear  stress  profiles  and  profiles  of  the  root  mean  square 
value  of  the  fluctuating  velocity  gradient  at  the  wave  surface.  The 
measurements  were  made  in  the  large  aspect  ratio  channel  used  by  Thorsness 
[70].  Wall  stresses  were  measured  with  electrochemical  techniques. 

To  test  the  theory  developed,  data  had  to  be  obtained  over  the  range 
of  dimensionless  wave  numbers  of  0.0005  <  a  =  <  0.01.  It  was 

possible  to  meet  the  above  requirement  by  modifying  existing  equipment 
in  one  of  two  ways.  The  facility  could  be  redesigned  so  as  to  increase 
the  range  of  friction  velocities  attainable  by  a  factor  of  four  or  waves 
of  four  times  the  previously  used  wavelength  could  be  employed.  Cook  [15] 
has  shown  that  for  a  given  channel  the  minimum  channel  height  to  wave¬ 
length  ratio  that  can  be  used  before  the  upper  wall  interferes  with  the 
measurements  at  the  wave  surf ace  is  £  -  1.  This  condition  would  have 
required  the  construction  of  a  new  channel  if  waves  of  longer  wavelength 
were  to  be  used.  For  this  reason  it  was  decided  to  use  the  existing 
channel  and  to  install  a  pumping  network  that  would  provide  the  necessary 
range  of  friction  velocities. 

Since  the  measurements  are  to  be  compared  with  linear  theory,  the 
amplitude  of  the  wave  model  used  in  gathering  the  data  had  to  be 
sufficiently  small  so  as  to  prevent  the  appearance  of  higher  order 

Zilker  has  demonstrated  that  provided  the 


harmonics  in  the  data. 
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dimensionless  amplitude, 


a  -  u 

CL 


<  30,  the  contribution  of  second  and  higher 


order  harmonics  to  the  total  shear  stress  profile  is  negligible. 


(a)  Flow  Loop 


The  flow  loop  used  in  this  study  was  originally  built  by  G.  W.  Cook 
and  modified  by  D.  P.  Zilker.  The  modifications  necessary  to  carry  out 
the  experiments  reported  in  this  thesis  are  reported  in  detail  by 
J.  J.  Buckles  [9].  Consequently  the  description  given  here  is  limited 
to  the  essential  features. 

The  flow  loop  shown  in  Figure  4.1  consists  of  a  rectangular  channel 
having  a  cross  sectional  area  of  2  in.  x  24  in.  and  a  length  of  27.5  ft. 

The  channel  was  originally  designed  to  handle  liquid  flows  for  which 
the  velocity  profiles  could  be  considered  fully  developed  in  the  neigh¬ 
borhood  of  the  test  section.  This  portion  of  the  channel  consists  of  a 
24  in.  x  27  in.  removable  section  located  at  the  far  end  of  the  27.5  ft. 
run.  It  is  here  that  the  wave  model  used  in  the  experiments  was  inserted. 

Two  pumps  connected  in  parallel  were  used  in  the  experiments.  The 
smaller  pump  is  a  Worthington  6  CNG84  centrifugal,  316  stainless  steel 
model.  It  was  driven  by  a  5  h.p.  motor.  This  pump  was  used  to  deliver 

flow  rates  of  up  to  800  (g.p.m.)  corresponding  to  a  channel  Reynolds  number, 

Ubh/2 

— - —  =  42,000, at  room  temperature.  The  larger  pump  constructed  from 

316  stainless  steel  is  a  Worthington  model  1050-D  centrifugal  pump  and 

is  driven  by  a  60  h.p.  motor.  This  pump  is  able  to  deliver  flow  rates 

U,  h/2 

of  up  to  2700  (g.p.m.)  corresponding  to  — ^ —  -  120,000.  The  range  of 

-2 

friction  velocities  obtained  with  the  two  pumps  is  2.93  x  10  ft/sec  < 

*  -1 

u  <  5.86  x  10  ft/sec.  The  pumping  network  is  shown  in  Figure  4.2. 

All  piping  and  fittings  were  either  6  in.  or  8  in.  Celanese  schedule 


4^ 

O 


Figure  4.1.  Schematic  Diagram  of  Flow  Loop 

1.  Down  stream  rectangular  to  round  diffuser 

2.  Removable  wave  surface 

3.  Test  section 

4 .  Channel 

5 .  Honeycomb 

6.  Up  stream  rectangular  to  round  diffuser 

7.  Annubar  flow  meter 

8.  Butterfly  throttling  valve 

9.  Removable  blanking  plate 

10.  Diaphragm  valve 

11.  Small  pump 

12.  By  pass  diaphragm  valve 

13.  Large  pump 

14.  Cooling  coils 

15.  Reservoir  tank 
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Figure  4.1.  Schematic  Diagram  of  Flow  Loop 
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Figure  4.2.  Photograph  of  Pumping  Network 
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80  P.V.C.  Provisions  were  available  for  flow  rate  and  temperature 
control.  An  additional  set  of  cooling  coils  had  to  be  installed  in  the 
recirculating  tank  to  maintain  temperature  control  when  running  with 
the  large  pump.  The  flow  rate  was  monitored  by  using  an  annubar 
connected  to  either  a  mercury  or  mirriam  oil  filled  manometer. 

(b)  Test  Section 

2  a 

A  wave  with  a  -  ratio  of  0.014  and  a  wavelength  of  2  in.  was 

A 

used  in  this  study.  Six  cutting  tools  were  made  for  the  fabrication  of 

the  wave  surface.  Four  of  the  cutting  tools  were  used  to  construct 

complete  waves.  Figure  4.3,  while  the  remaining  two.  Figure  4.4,  were 

used  so  as  to  mesh  the  leading  and  lagging  waves  into  the  flat  portion 

of  the  wave  section  which  lined  up  with  the  channel.  The  cutting  tools 

were  made  from  hardened  steel  using  a  new  approach.  A  computer  controlled 

mill  which  automatically  compensated  for  the  radius  of  curvature  of  the 

mill  was  used  to  generate  the  tools.  The  tools  were  2  in.  in  length 

and  machined  according  to  the  relationship  y  =  a^  sin  (2t t  x/A).  The 

cutting  edge  of  the  tools  had  a  9*5°  slope  and  therefore  the  amplitude 

of  the  tools  had  to  be  compensated  for  in  order  to  provide  the  required 

wave  amplitude.  The  formula  used  to  calculate  this  compensation  is 

a, (actual  wave)  =  a,  (cutting  tool) /cos  (8),  where  9  is  the  angle  of  the 
d  d 

cutting  tool  edge. 

The  tools  were  placed  in  a  rotary  cutter  assembly  and  a  wave  train 
consisting  of  ten  complete  wavelengths  was  machined  into  a  27  in.  x  24  in.  x 
2  in.  thick  plexiglass  section.  The  waves  were  machined  perpendicular  to 
the  mean  flow  direction.  Finally  the  two  waves  which  meshed  into  the 
flat  portions  of  the  wave  section  were  machined.  A  cross  sectional  view 
of  the  wave  pattern  is  shown  in  Figure  4.5. 
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Amplitude  of  Cutting  Tool  Exaggerated 


9-1/2° 


~H  H-l/4" 


Figure  4.3.  Cross  Sectional  View  of  Cutting  Tool 


Amplitude  of  Tool  Exaggerated 


Figure  4.4.  Cross  Sectional  View  of  End  Cutting  Tool 


Amplitude  of  Waves  Exaggerated 
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Figure  4.5.  Cross  Sectional  View  of  Wave  Surface 
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Holes  for  the  electrodes  were  drilled  into  the  wave  surface 
according  to  the  following  pattern.  Twenty-one  holes  of  0.025  in. 
diameter  were  drilled  into  the  third  and  sixth  waves  respectively  while 
a  combined  total  of  forty-one  holes  were  drilled  into  the  eighth  and 
ninth  waves.  The  electrode  pattern  shown  in  Figure  4.6  was  chosen  so 
as  to  minimize  the  interference  between  the  electrodes  and  to  test 
whether  a  fully  developed  flow  field  existed. 

Platinum  electrodes,  0.020  in.  diameter,  were  epoxied  into  each 
of  the  holes  using  Techkits  A-12  epoxy  as  follows.  The  epoxy  was 
injected  into  a  hole.  Before  the  epoxy  dried,  a  two  inch  piece  of 
platinum  wire,  soldered  onto  a  twelve  inch  coated  copper  wire,  was 
inserted  into  the  hole  from  the  back  of  the  wave  section.  A  small 
amount  of  platinum  wire  was  allowed  to  protrude  on  the  wave  side.  The 
platinum  wire  was  coated  with  a  thin  layer  of  epoxy  prior  to  insertion 
into  the  holes. 

Once  dry  the  electrodes  were  filed  down  flush  with  the  wave  surface 
using  a  soft  file.  This  had  the  advantage  of  not  damaging  the  plexi¬ 
glass  while  removing  the  excess  platinum.  Once  the  electrodes  were 
flush,  the  surface  was  sanded  down  with  progressively  finer  grades  of 
sandpaper,  polished  with  DUPONT  0861N  Rubbing  Compound,  and  DUPONT 
0761N  Polishing  Compound,  and  finally  polished  with  Mirror  Glaze  Plastic 
Cleaner  and  Polish.  On  completion  of  the  polishing^  the  profile  of  the 
waves  was  checked  using  a  dial  indicator. 

The  shear  stress  measurements  were  made  over  the  sixth,  eighth  and 
ninth  waves.  Figure  4.7,  4.8  and  4.9  show  the  profiles  of  the  above 
waves  respectively  as  a  function  of  x/X,  taking  the  wave  trough  to  be 
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Figure  4.6.  Electrode  Pattern 
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Figure  4.7.  Profile  of  Wave  6 
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Figure  4,8.  Profile  of  Wave  8 
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Figure  4.9.  Profile  of  Wave  9 
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zero.  The  amplitude  of  each  wave  was  estimated  by  performing  a  least 
squares  analysis  of  the  wave  profile  data.  The  maximum  relative  phase 
shift  among  the  three  wave  profiles  is  less  than  3  .  The  tabulated  wave 
profile  data  is  listed  in  Appendix  F.  Figures  4.10  and  4.11  show  the 
front  and  back  view  of  the  completed  wave  section. 


(c)  Electrochemical  Technique 

The  electrochemical  technique  which  uses  a  diffusion  controlled 
electrode  for  the  measurement  of  wall  shear  stress  was  developed  by 
Reiss  [58]  and  Mitchell  and  Hanratty  [47].  A  chemical  reaction  occurs 
on  an  electrode  which  is  embedded  flush  with  the  surface.  The  resulting 
reaction  current  is  related  to  the  wall  shear  stress.  The  electro¬ 
chemical  reaction  employed  is  an  oxidation  reduction  couple  employing 
a  large  excess  of  supporting  electrolyte.  The  cathode  which  is  the 
shear  stress  probe  is  much  smaller  than  the  anode  which  makes  the  cathodic 
reaction  the  limiting  step.  By  suitably  adjusting  the  applied  potential 
the  concentration  of  active  specie  at  the  surface  can  be  made  equal  to 
zero.  The  result  is  that  the  reaction  rate  and  hence  current  flow 
becomes  a  function  only  of  the  rate  of  diffusion  to  the  cathode.  Con¬ 
sequently  the  current  flow  can  be  related  directly  to  the  mass  flux  at 
the  surface  and  the  physical  properties  of  the  system. 

The  redox  reaction  used  in  this  study  is  the  potassium  iodide  and 
iodine  reaction  system,  in  which  the  following  reactions  occur. 

l“  +  2e~  - *■  3I~  (cathode) 


31* 


I~  +  2e 


(anode) 


The  approximate  concentrations  of  I3  and  the  potassium  iodide  supporting 
electrolyte  were  0.0015m  and  0.2m  respectively.  The  properties  of  the 
electrolyte  are  summarized  in  Appendix  G. 


Figure  4.10.  Photograph  of  Front  View  of  Wave  Section 
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Under  these  conditions  a  mass  transfer  coefficient  may  be  defined 


by  the  relationship 
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FAC, 
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(4.1) 


where  I  is  the  current  at  the  cathode, ne  is  the  number  of  electrons 
involved  in  the  reaction,  F^  is  Farraday's  constant,  Ag  is  the  area  of 
the  electrode  and  C,  is  the  bulk  electrolyte  concentration.  Through  a 

D 

solution  of  the  mass  balance  equations  Mitchell  and  Hanratty  [47]  have 
shown  that  the  instantaneous  mass  transfer  coefficient  is  related  to 
the  wall  shear  stress  as  follows 


f  3  1 

2T(4/3) K 
m 

»v 

3 

w  > 

D2 

j 

(4.2) 


where  y  is  the  viscosity  of  the  fluid,  D  is  the  diffusion  coefficient  for 
the  reacting  species  andLgis  the  equivalent  length  of  the  electrode, 
equal  to  0.816  times  the  diameter  of  the  electrode.  Using  equations 
(4.1)  and  (4.2)  and  neglecting  the  transverse  component  of  the  fluctu¬ 
ating  stress  x '  the  current  measured  can  be  related  to  the  wall  shear 
zd 

stress  as  follows 


i  =  c  ( xw  +  t;  ) 


1/3 


X  , 


(4.3) 


where  C  is  a  proportionality  constant  and  tw  and  x^  are  the  time  averaged 

d  d 

and  fluctuating  wall  shear  stresses  respectively.  Furthermore  if 
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(4.5) 


and 


I  -  I 

T 


(4.6) 


Equations  (4.5)  and  (4.6)  are  the  design  equations  used  to  determine 
the  average  and  root-mean  square  value  of  the  wall  shear  stress  respectively. 

(d)  Data  Acquisition 

The  electronic  equipment  used  to  measure  the  electrode  currents  was 
the  same  as  that  described  in  previous  work,  (Zilker  [77]).  Figure  4.12  gives 
a  schematic  representation  of  the  basic  electronics.  A  118  A  amplifier 
was  used  to  apply  a  negative  potential  to  the  alternate  electrodes  and 
to  act  as  a  current  to  voltage  converter.  The  applied  potential  was 
set  by  the  adjustment  of  a  Helipot  so  that  the  potential  was  always  kept 
in  the  mass  transfer  controlled  plateau  region.  A  feedback  resistor, 

R^,  was  used  such  that  the  output  voltage,  V  was  between  2.5  and  3.5  volts. 
The  output  voltage  was  sampled  with  an  Isaac-Cyborg  twelve  bit  analog  to 
digital,  A/D,  converter.  The  A/D  converter  was  in  turn  linked  to  an 
Apple  II  plus  computer.  Up  to  sixteen  channels  could  be  sampled  simul¬ 
taneously  at  a  sampling  rate  of  1  KH^.  A  block  diagram  representation 
of  the  data  acquisition  process  is  shown  in  Figure  4.13. 

Data  was  stored  in  an  integer  basic  format  which  allowed  the  Apple  II 
to  store  approximately  22,000  data  points  in  memory.  A  6502  assembly 
language  program  was  written  for  the  purpose  of  analyzing  the  data.  The 
program  is  able  to  calculate  the  mean  and  the  variance  of  the  output 
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Figure  4.12.  Electronic  Circuit  for  Shear  Stress  Measurements 


Computer 


Figure  4.13.  Schematic  Diagram  of  Data  Acquisition  System 
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voltage  minus  the  applied  voltage, (Vq  -  Vapp) .  Using  these  values  and 
the  design  equations,  (4-5)  and  (4.6)  the  average  shear  stress  and  the 
root  mean  square  value  of  the  fluctuating  shear  stress  can  be  obtained 
as  described  in  Appendix  H.  Provision  was  also  made  to  transfer  the  raw 
data  to  a  floppy  disc  so  as  to  allow  for  further  analysis  if  necessary. 

(e)  Shear  Stress  Measurements 

Shear  stress  measurements  were  obtained  using  the  previously  described 
wave  surface.  Before  each  run  was  made  the  wave  surface  was  cleaned  with 
soapy  water  using  a  soft  sponge.  This  was  followed  by  thoroughly  rinsing 
the  wave  surface  with  deionized  water.  The  channel  was  filled  with 
electrolyte  and  allowed  to  run  for  a  period  of  time  in  order  to  remove 
any  entrapped  bubbles.  The  temperature  was  maintained  at  (26  ±  0.2)  C  . 

At  the  conclusion  of  a  run  a  sample  of  fluid  was  taken  and  the  required 
chemical  and  physical  tests  were  performed  as  outlined  in  Appendix  G. 

Figures  4.14  through  4.34  show  the  measured  time  average  shear  stress 
profiles  of  selected  number  of  runs.  Measurements  of  the  root  mean  square 
level  of  the  fluctuating  shear  stress  are  shown  in  Figures  4.35  through 
4.40.  The  average  and  fluctuating  shear  stress  measurements  are 
normalized  with  the  wave  averaged  shear  stress.  The  data  are  plotted 
with  respect  to  a  cosine  wave  such  that  x/X  =  0  corresponds  to  the  upstream 
crest  and  x/X  =  1  corresponds  to  the  downstream  crest.  On  each  plot  a 
line  is  drawn  through  the  points  representing  the  best  sinusoidal  fit 
of  the  data.  The  curves  were  obtained  by  performing  a  least  squares 
analysis  of  the  data.  The  details  of  this  analysis  are  presented  in 
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Figure  4.14.  Shear  Stress  Distribution  for  Re  =  5,970 
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Figure  4.16.  Shear  Stress  Distribution  for  Re  =  8,450 
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Figure  4.17. 


Shear  Stress  Distribution  for  Re  =  9,650 
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Figure  4.18.  Shear  Stress  Distribution  for  Re  =  13,000 
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Figure  4.19*  Shear  Stress  Distribution  for  Re  =  14,600 
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Figure  4.20.  Shear  Stress  Distribution  for  Re  =  15,700 
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Figure  4.21. 


Shear  Stress  Distribution  for  Re  =  19,800 
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Figure  4.22.  Shear  Stress  Distribution  for  Re  =  20,600 
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Figure  4.23.  Shear  Stress  Distribution  for  Re  =  27,950 
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Figure  4,24.  Shear  Stress  Distribution  for  Re  =  30,620 
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Figure  4.25.  Shear  Stress  Distribution  Re  =  34,700 
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Figure  4.26.  Shear  Stress  Distribution  for  Re  -  36,000 
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Figure  4.27.  Shear  Stress  Distribution  for  Re  =  58,000 
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Figure  4.28.  Shear  Stress  Distribution  for  Re  =  64,000 


u> 


Wove  Height 

(in.) 


X 


Figure  4.29.  Shear  Stress  Distribution  for  Re  «  71,000 
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Figure  4.30.  Shear  Stress  Distribution  for  Re  =  73,500 
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Figure  4.31.  Shear  Stress  Distribution  for  Re  «  84,500 
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Figure  4.32.  Shear  Stress  Distribution  for  Re  =  102,500 
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Figure  4.33.  Shear  Stress  Distribution  for  Re  =  114,500 
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Figure  4.34.  Shear  Stress  Distribution -for  Re  =  122,500 
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Figure  4.35.  Fluctuating  Shear  Distribution  for  Re  =  6,680 
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Figure  4,36.  Fluctuating  Shear  Stress  Distribution  for  Re  =  8,450 
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Figure  4.37.  Fluctuating  Shear  Stress  Distribution  for  Re  =  9,650 
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Figure  4.38.  Fluctuating  Shear  Stress  Distribution  for  Re  =  13,000 
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Figure  4,39.  Fluctuating  Shear  Stress  Distribution  for  Re  =  15,700 
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Figure  4,40.  Fluctuating  Shear  Stress  Distribution  for  Re  =  19,800 
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Figures  4.41  and  4.42  summarize  the  average  shear  stress  data  by 
plotting  the  phase  angles  and  amplitudes  of  the  sinusoidal  curves 
fitted  to  the  wall  shear  stress  measurements  as  a  function  of  Reynolds 
number.  The  phase  angle  is  measured  with  respect  to  the  wave  surface, 
with  shifts  in  the  direction  opposite  the  flow  direction  taken  as 
positive.  The  magnitude  of  the  variation  of  shear  stress  decreases  with 
Reynolds  number  while  the  phase  angle  increases  initially,  passes  through 
a  maximum  at  about  a  Reynolds  number  of  35,000  and  then  decreases. 

Figures  4.43  and  4.44  summarize  the  phase  angles  and  amplitudes  of 
the  fluctuating  shear  stress  measurements.  An  error  in  the  amplitudes 
of  these  measuremtns  is  expected  due  to  the  spatial  averaging  of  the 
velocity  fluctuations  over  the  surface  of  the  electrode  and  due  to  the 
problem  of  frequency  response.  Mitchell  and  Hanratty  [47]  have  shown 
how  to  correct  for  both  of  these  effects.  Figure  4.47  also  shows  the 
corrected  amplitudes.  The  details  of  the  corrections  are  given  in 
Appendix  J.  The  magnitudes  of  the  fluctuations  decrease  with  increasing 
Reynolds  number  while  the  phase  angles  increase  with  Reynolds  number. 

The  shear  stress  and  intensity  data  are  tabulated  in  Appendix  F. 


Figure  4.41,  Phase  Angle  of  Average  Shear  Stress  Measurements  as 
a  Function  of  Reynolds  Number 
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Figure  4.42.  Amplitude  of  Average  Shear  Stress  Measurements 
a  Function  of  Reynolds  Number 
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Figure  4.43.  Phase  Angle  of  Fluctuating  Shear  Stress  Measurements  as  a  Function  of  Reynolds  Number 
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Figure  4.44.  Amplitude  of  Fluctuating  Shear  Stress  Measurements  as  a 
Function  of  Reynolds  Number 
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number  of  points  in  the  range  x^,  .  ..,  x  .  This  leads  to  a  system  of, 
n,  algebraic  equations  which  can  then  be  solved  for  the  coefficients,^. 
The  system  of  algebric  equations  will  be  linear  provided  the  differential 
equation  being  solved  is  linear.  Villadsen  and  Michelsen  [74]  provide 
a  detailed  review  of  this  method. 

Galerkin's  method  (Strang  and  Fix  [68])  is  based  on  the  orthogonality 
of  functions.  If  y(x)  represents  an  exact  solution  to 


dnf(x) 


dx 


n 


g(f(x) . fn(x),  v(x))  , 


(5.3) 


then  y(x)  has  an  associated  residual,  r(x),  defined  as 


r (x)  ■  d  -  g(y(x),  yn(x),  v(x>)  , 


(5.4) 


dx 


which  is  identically  equal  to  zero  for  all  x,  (since  y(x)  represents 
an  exact  solution).  Therefore,  the  residual,  r(x),  is  orthogonal  to 
every  function  and  in  particular  it  would  be  orthogonal  to  the  set  of 
basis  functions.  Hence 

•b 

r (x)  4>i(x)  dx  =  0  i  =  1,  .  .  .n.  (5.5) 

Ja 


Since,  y(x),  is  a  linear  combination  of  basis  functions  it  is  not 
expected  to  be  an  exact  solution  of  the  differential  equation  being 
solved.  What  GalerkinTs  method  does  is  to  choose  a  y(x)  which  has  a 
residual  that  is  orthogonal  to  all  basis  functions,  ^(x),  ...,  <J>^(x), 
respectively.  This  procedure,  which  requires  the  vanishing  of  a  number 
of  integrals,  again  leads  to  the  solution  of  a  system  of  algebraic 
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equations  for  the  unknown  constants,  c_.  .  Further  details  regarding 
the  above  projection  methods  can  be  found  in  Prenter  [55],  Lucas  and 
Reddisen  [43],  Russel  [6],  and  Ortega  and  Poole  [52]. 

The  next  class  of  methods  are  the  so  called  shooting  methods  where 
the  basic  aim  is  to  always  solve  an  initial  value  problem  by  using  any  one 
of  the  standard  forward  integration  techniques.  This  is  achieved  by 
guessing  any  parameters  which  are  needed  to  make  this  possible,  and  by 
adjusting  the  guesses,  directly  with  linear  problem  or  iteratively  with 
nonlinear  problems,  such  that  all  the  boundary  conditions  are  satisfied. 
There  are  many  variations  of  the  basic  idea.  These  include  direct 
shooting  or  the  method  of  superposition  for  linear  problems,  (Scott  and 
Watts  [  63],  Na  [  50]).  For  nonlinear  problems  direct  shooting  must  be 
used  in  conjunction  with  Newtons  method,  (Na  [50]).  Alternatively 
several  workers  including  Bellman  and  Kalaba  [  5  ]  and  Scott  and  Watts 
[62  ]  have  suggested  using  the  method  of  superposition  together  with 
quasilinearization. 

In  the  final  class  of  methods  low-order  finite  difference  formulae 
are  used  and  applied  as  approximations  to  the  differential  equation  at 
a  number  of  discreet  points  in  the  range.  The  boundary  conditions  are 
satisfied  exactly  if  they  do  not  involve  derivatives  or  approximately 
if  they  do.  The  result  of  approximating  the  differential  equation  by 
finite  difference  formulae,  is  a  set  of  linear  algebraic  equations  in 
the  case  of  a  linear  differential  equation  and  a  set  of  nonlinear  algebraic 
equations  in  the  case  of  a  nonlinear  differential  equation.  The  method  is 
then,  to  solve  the  resulting  set  of  algebraic  equations  directly  for 
linear  problems  and  iteratively  for  nonlinear  problems,  to  give 
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approximations  to  the  solution  of  the  differential  equation  simultaneously 
at  all  mesh  points.  The  method  had  been  reviewed  by  Ortega  and  Poole 
[52]  and  Fox  [ 20  3  . 

The  linear  differential  equations  (3 . 13)- (3 . 17)  and  (3. 72)  —  (3. 77) 
which  specify  the  wave  induced  stream  function,  F(y),  for  the  Zero  Equation 
and  K-e  Models  were  solved  using  a  modified  version  of  the  method  of  super¬ 
position  to  be  discussed  in  the  next  section.  The  main  advantage  of  the 
method  of  superposition  (direct  shooting)  is  the  existence  of  very  good 
and  well  programmed  methods  for  solving  systems  of  initial  value  problems. 
These  integrators  are  programmed  to  adjust  the  step  by  step  intervals  so 
as  to  provide  the  desired  level  of  accuracy  required  by  the  user. 

The  nonlinear  differential  equations  (3 . 67)- (3. 71)  which  specify 
the  mean  flow  field  when  using  the  K-e  Model  were  solved  by  combining 
finite  difference  techniques  with  quasilinearization,  to  be  discussed  in 
a  subsequent  section.  Initial  attempts  were  made  to  solve  this  problem 
using  direct  shooting  and  Newton’s  method.  However,  due  to  the  nature 
of  the  equations,  if  the  initial  conditions  are  not  known  accurately  the 
solution  blows  up  very  fast.  An  alternate  method  suggested  by  Scott  [64], 
but  not  tried,  is  the  method  of  parallel  shooting. 

(a)  Numerical  Solution  of  the  Wave  Induced  Flow 

A  linear  boundary  value  problem  such  as  equations  (3 . 13)- (3. 17)  or 
equations  (3 . 72)- (3 . 77)  may  be  written  as  a  system  of  first  order  equations 
and  expressed  in  the  following  form,  (see  Appendix  K) , 


y’  (x)  =  F(x)  y (x)  +  g(x) , 


(5.6) 
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Ay  (a)  =  a  , 

(5.7) 

By(b)  =  3  , 

(5  8) 

where  y(x)  is  a  column  vector  of  length  n,  F(x)  is  an  n  x  n  matrix  of 
rank  k,  A  is  an  (n  -  k)  x  n  matrix  of  rank  (n  -  k)  ,  B  is  k  x  n  matrix 
of  rank  k,  and  a  and  6  are  column  vectors  of  (n  -  k)  and  k  components 
respectively. 

Applying  the  method  of  superposition  any  solution  of  equations 
(5.6)- (5. 8)  can  be  expressed  as 

y (x)  =  v(x)  +  Cjy^x)  +  c2y2(x)...  ckyk(x)  (5.9) 

=  v  (x)  +  U  (x)c  (5.10)- 

where  c  is  a  vector  with  k  components,  v(x)  is  a  solution  of  the 
nonhomogeneous  system 

v * (x)  =  F(x)  v (x)  +  g(x),  (5.11) 

Av(a)  *  a  (5.12) 

and  U(x)  is  an  n  x  k  matrix  whose  columns  yJ_(x)...,  y  (x)  are  solutions 
of  the  homogeneous  system 


U!  (x)  =  F(x)  U(x)  ,  (5.13) 

AU(a)  =  0.  (5.14) 


The  constant  vector  c  is  determined  by  satisfying  the  final  boundary 
condition  at  x  =  b,  equation  (5-8). 


By (b)  -  BU (b) c  +  Bv(b)  =  @ 


(5.15) 
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It  should  be  mentioned  that  in  classical  superposition  n  homogeneous 
solutions  would  have  to  be  generated.  However  by  suitably  choosing  the 
initial  conditions,  equations  (5.12)  and  (5.14),  it  is  possible  to 
generate  only  k  homogeneous  solutions.  If  the  direction  of  integration 
is  reversed  (i.e.  from  b  to  a)  then  (n  -  k)  solutions  are  required.  For 
simple  boundary  conditions  [equations  (3.16),  (3.17)  and  (3.77)]  the 
initial  matrix,  U(a),  can  be  chosen  by  inspection  (see  Appendix  K) .  For 
complex  boundary  conditions  Scott  and  Watts  [62]  have  developed  a 
general  procedure  for  generating  the  initial  matrix,  U(a). 

It  often  happens  that  the  exact  mathematical  procedure  for  obtaining 
the  solution  of  equations  (5. 6)- (5. 8)  leads  to  poor  or  even  completely 
incorrect  results  when  applied  as  a  numerical  procedure.  This  occurs 
when  the  matrix  F(x)  has  eigenvalues  whose  real  parts  are  well  separated. 
In  this  case  the  various  homogeneous  solutions  and  the  particular  solution 
all  grow  at  vastly  different  rates.  As  a  consequence  of  the  finite  word 
length  associated  with  digital  computers  the  rapidly  growing  solutions 
eventually  swamp  all  other  solutions,  making  them  proportional  to  these 
growing  solutions  or  combinations  of  them.  In  this  case  the  matrix  U  will 
become  increasingly  ill  conditioned  such  that  the  vector  c  determined  by 
equation  (5.15)  will  yield  incorrect  results.  The  above  mentioned 
problem  arises  in  the  solution  of  equations  (3. 13)-(3. 17)  and  (3.72)- 
(3.77). 

Consider  equations  (3.13)  and  (3*14)  at  large  y  such  that  the 
derivatives  of  the  mean  velocity  profile  Ufand  U"  are  small  and  the 
turbulence  term,  R  ,  is  negligible.  Then 

F""  -  (2a2  +  iaU)  F"  +  (a4  +  ia3  U)  F  -  ia3U  2  -  a4U  =  0 


(5.16a) 
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Upon  application  of  the  standard  method  of  solution  for  linear 
ordinary  differential  equations  with  constant  coefficients  and  keeping 
the  two  solutions  which  are  bounded  at  large  y  there  results 

p  -  \ + \ + fp  ■  «-i6b> 


where 


exp  (-ay)  +  A 2 

_  1/2 
+  iaU) 


exp  (-By) 


+  U, 


(5.16c) 

(5.16c) 


Assuming  infinite  precision  the  two  homogeneous  solutions  exp  (-ay) 
and  exp  (-By)  could  be  generated  by  solving  the  homogeneous  form  of  (5.16a) 
subject  to  the  initial  conditions 

fH;l  =  1  »  *2  =  (-“)n  y  =  (5.16d) 

X  *  1>  ^2  =  ("6)n  7  =  ^  (5 . 16e) 

where  Fn  denotes  the  nth  derivative. 

In  fact  these  initial  conditions  will  produce  the  following 
homogeneous  solutions: 


Fr  =  exp  (ay^)  exp  (-ay) 


(5.17a) 


Fr  =  exp  (By^)  exp  (-By)  . 


(5.17b) 


Since  computers  do  not  have  infinite  precision  let  the  total  amount 
of  round  off  and  truncation  error,  e,  (Lightfoot  [41])  enter  the  problem 
in  the  initial  conditions  as  follows: 
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„  n  ,  x n 

F„  =  1  +  e  F  =  (-a)  at  yo 

rtl  1 


_  -  1  +  e  F"  =  (-8) 

FH2  H2 


at  y 

J  nr 


(5.18a) 

(5.18b) 


The  numerical  solutions  that  would  be  generated  using  the  above 
boundary  conditions  are: 


H, 


Be 


ct  Be 


1  + 


2(8  -  a)  2(82  -  a2) 


eay~  e‘ay 


a  B  e 


a  Be 


2a(a  +  B)  (B2  -  a2) 


-ay  ay 
e  e 


eey„  e-ey  . — 11—  e“ey~  eBy  ,  (5.i9a) 


2  2 

2(fT  -  a  ) 


2  2 
2(B  -  ct) 


Hn 


ae 


a  B  e 


1  + 


2(o-6)  2(o2-62) 


eSy"  e-6y 


a  8  e  a  8  e 

-  + 


26(0  +  6)  (o2  -  62) 


e-sy»  e8y 


2  2 
Be  Be 

_ _  eayoo  e“«y  _  z _  e’ayco  e«y 


2  2 
2(aZ  -  B  ) 


2  2 
2(aZ  -  B  ) 


(5.19b) 


The  behavior  of  F„  and  F  as  y — >0  yields 
H1  H2 
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-g2e  Jy=o  -By 

2  2’ 
2 (gZ  -  a) 


(5.20a) 


F  eeyoo  e-By  #  (5.20b) 

H2 

since  |  3  |  >>  |  a  |  . 

Hence  F  and  F  are  no  longer  independent  solutions  of  the  homo- 
H1  H2 

geneous  form  of  equation  (5.16a).  The  point  at  which  loss  of  independence 
occurs  depends  on  the  word  size  of  the  computer  being  used. 

Several  workers  have  recognized  this  difficulty  and  have  proposed 
methods  for  maintaining  the  linear  independence  of  the  solutions.  These 
include  Kaplan's  filtering  technique  [22],  Davey's  [17]  complete  ortho¬ 
normalization  procedure,  and  the  method  of  compound  matrices  Ng  and 
Reid  [51] .  The  most  widely  used  and  tested  method  is  the  Gram-Schmidt 
orthonormalization  procedure. 

The  use  of  the  Gram-Schmidt  orthonormalization  procedure  was  first 
proposed  by  Godunov  [25]  and  independently  by  Bellman  and  Kalaba  [5]. 

The  method  has  been  discussed  and  used  by  Conte  [14],  Gersting  [23] 
and  Scott  and  Watts  [63].  Scott  and  Watts  [62]  have  refined  the  basic 
method  and  incorporated  these  modifications  into  an  excellent  code 
named  Suport  [62]. 

In  this  thesis  the  Suport  code  was  used  together  with  an  indepen¬ 
dently  written  code.  The  latter  was  developed  to  obtain  a  better 
understanding  of  the  orthonormalization  procedure.  The  advantage  of 
the  Suport  code  is  its  extreme  efficiency  and  the  modifications  which 
have  been  added  that  take  advantage  of  the  complex  structure  of  problems 
like  the  Orr — Sommerfeld  equation  (see  Watts,  Scott  and  Lord  [75]). 
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(b)  Orthonormalization  Procedure 

The  orthonormalization  procedure  outlined  here  is  the  method 
suggested  by  Scott  and  Watts  [63  ],  Conte  [14  ]  and  Godunov  [25  ]. 

At  the  points  where  orthonormalization  is  necessary  the  previous 
independent  set  of  base  solutions  are  converted  into  a  new  orthonormal 
set.  The  orthonormalization  is  achieved  by  employing  the  "modified" 
Gram-Schmidt  precedure  described  below.  Scott  and  Watts  [62  ]  have 
suggested  that  the  modified  Gram-Schmidt  procedure  be  used  since  it  is 
(1)  more  stable  (2)  easier  to  code  (3)  more  economical  on  storage. 
Thus  at  any  point  where  orthonomalization  is  judged  necessary 


U  -  PU 
old  new 


(5.21) 


where  U  ,,  is  the  old  independent  set,  U  is  the  new  orthonormal 
old  new 

basis  and  P  is  an  upper  triangular  matrix  (defined  in  the  next  section) 

that  is  constructed  from  the  columns  of  U  -  , .  At  each  orthonormalization 

old 

point  the  particular  solution  is  turned  into  the  orthogonal  complement 
of  the  base  set 


v 

new 


Vold 


-  U  w 
new  new 


(5.22) 


where  w  is  defined  in  the  next  section.  The  terms  v  , ,  and  v 

new  old  new 

refer  to  the  old  and  new  particular  solutions  respectively. 

Let  the  points  at  which  orthonormalization  has  been  carried  out  be 

denoted  by  x9,  x  and  let  x  =  a.  Now  consider  the  first 

-L  n  o 

interval  defined  in  the  range  [xq,  x  ].  Let  the  particular  solution, 

homogeneous  solutions  and  overall  solution  be  denoted  by  v  (x) .  U  (x) 

o  ’  o 

and  yQ(x)  repsectively.  Then  vq(x)  satisfies 
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v 1 (x)  =  F (x)v  (x)  +  g(x),  Av(a)  =  a  ,  (5.23) 

o  o 

and  U  (x)  satisfies 
o 

lT(x)  =  F(x)U  (x),  AUo(a)  =  0  .  (5.24) 

The  general  solution  on  this  interval  is  obtained  using  the  method  of 
superposition.  Hence, 

y  (x)  =  v  (x)  +  U  (x)  c  (5.25) 

J  o  o  o  o 

where  c  is  a  constant  to  be  determined  by  matching  y  (x  )  with  y  (x  ) . 

O  O  X  i.  x 

In  using  this  representation  Uq(x^)  and  vq(x^)  are  t^ie  °ld  independent 
set  of  base  vectors  and  the  old  particular  solution  at  x^.  The  new 
orthonormal  set  of  base  vectors  and  orthogonal  complement  to  this  set 
at  x^  are  denoted  by  u1(x1)  and  v^x^)  respectively.  Using  U^xp  and 
vl(Xi)  as  initial  conditions  on  the  interval  (x^,  x^)  the  above  procedure 
is  repeated  by  solving 

v^(x)  =  F(x)v1(x)  +  g(x)  (5.26) 

subject  to  the  initial  condition 

vl(Xi)  -  vq(x]l)  -  U1(x1)w1  (5.27) 

and  by  solving  the  system 

U^(x)  =  F (x)U^ (x)  (5.28) 


subject  to  the  initial  condition 
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W  ■  VX1)P1 


(5.29) 


The  general  solution  on  this  interval  is 


y^x)  =  v1(x)  +  U1(x)c1 


(5.30) 


where  the  vector  is  determined  by  matching  y^(x^)  with  y2(x2). 


yi(x2)  =  y2(x2) 


(5.31) 


Hence,  on  the  interval  (x^  x^+^)  the  solution  is  bbtained  as 


yi^Xi^  =  vi^  +  Ui^ci 


(5.32) 


provided  that  the  condition  of  continuity  is  met 


WV  ‘  yi(xi)  • 


(5.33) 


Thus  from  the  above  condition  of  continuity 


Vi_i(Xi)  +  =  vi^xi^  +  Ui^xi^ci  (5.34) 


but 


vi(Xi)  =  -  Ui(x±)wi 


(5.35) 


and 


“i-i^i)  ■  ui<xi>pi  • 


(5.36) 
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Therefore  equation  (5.34)  becomes 


vi-i(Ki>  +  ui<xi,pici-i 


+  0.(x.)(cl  -  w±) 


(5.37) 


This  provides  the  basic  recursion  formula  for  the  vector  c^_^  in 

terms  of  P . ,  w,  and  c. 

i  i  l 

Pici_1  =  (c±  -  w±)  (5.38) 

At  the  final  point,  x  =  b,  on  the  interval  (x^,  x^)  the  solution 
of  the  differential  equation  must  satisfy 


y 


m 


(b)  =  v  (b)  +  U  (b)c  . 
m  mm 


Hence,  using  equation  (5.8)  the  constant  c^  can  be  determined, 

By  (b)  =  Bv  (b)  +  BU  (b)c  =  8 
m  m  mm 


(5.39) 


(5.40) 


or 


BU  (b)c  =  8  -  Bv  (b)  .  (5.41) 

mm  m 

Having  determined  the  constant  vector  c^  (and  providing  all  the 

orthonormalization  information,  such  as  P^,  P^,  . Pm  have  been 

retained)  the  vectors  c  c  0,  . c  can  be  found  by  solving  the 
7  m-1  m-2*  o 

systems,  equation  (5.38), 

Pc  ^  =  c  -  w  r  =  m,  m-1,  ...,  1  (5.38) 

r  r-1  r  r 

by  back  substitution.  If  the  initial  values  IKCx^,  v^x^)  i  =  0,  1, 

2,  ...,  m  have  also  been  retained  the  solution  y(x)  can  be  calculated 
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at  any  point  xe  (x^,  xi+j_)  by  integrating  forward  from  and  using 
the  relation 


y(x)  =  vi(x)  +  Ui(x)ci  . 


(5.42) 


In  the  formulation  presented  above  the  particular  solution  is  not 
normalized  (i.e.  not  made  into  a  unit  vector).  Scott  and  Watts  [62] 
have  shown  that  the  final  solution  representation  with  and  without 
normalizing  the  particular  solution  are  mathematically  equivalent.  If 
the  particular  solution  is  normalized  then  the  recursion  formula  equation 
(5.38)  is  not  valid.  Scott  and  Watts  [62]  have  developed  a  recursion 
formula  valid  for  this  case.  They  mention  that  a  possible  advantage  of 
normalizing  the  particular  solution  is  that  it  is  then  well  scaled. 


(c)  Gram-Schmidt  Orthogonalization  Procedure 

In  this  section  the  modified  Gram-Schmidt  procedure  is  outlined 
together  with  a  set  of  criteria  that  must  be  met  to  prevent  the  solution 
vectors  from  becoming  dependent. 

Scott  and  Watts  [62]  have  suggested  that  the  old  independent  set 
of  vectors  Y  ^  be  reorganized  such  that  the  vector  with  the  largest 
norm  appears  in  the  first  column.  This  involves  postmultiplying  the 
vector  Y  by  a  pivot  matrix  E  such  that 


Y  .  ,  =  Y! 1  ,  E  -  U  P 
old  old  new 


(5.43) 


where  Y  is  the  reorganized  matrix.  The  term  Y^^  is  the  matrix 
obtained  directly  from  integration,  E  is  the  permutation  matrix,  U 
is  the  orthonormalized  basis  and  P  is  an  upper  triangular  matrix  defining 
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the  transformation.  If  the  permutation  matrix  is  included  in  the 
algorithm  then  the  recursion  formula  equation  (5.38)  must  be  modified. 
Consider  the  interval  (x.j_^,  x^ . 


yi-i(xi)  =  vi-i(xi} 


+  U±_1(xi)  E._l2i-1 


(5.44) 


where 


■  di(xi)p 


(5.45) 


Also 


yi(xi)  =  vi(xi)  +  ui(xi)  Ei£i  5 


-  vi(xi)  +  U±(xi)c1  . 


(5.46) 

(5.47) 


Equating  (5.44)  and  (5.47)  at,  x^  yields 


(5.48) 


Eliminating  v±  (x^  and  v±(xi),  the  following  recursion  relation  is 
obtained : 


Pi  £i-l  =  (°i  ~  Wi)  ’ 


(5.49) 


where 


c .  -  =  E.  .  c  - 
i-l  l-l  i-1 


(5.50) 


The  modified  Gram-Schmidt  procedure  defines  the  matrix  P  (see 


Stewart  [67])  as. 
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P  = 


(0) 


(0) 

U2  ’  U1 


.(1) 


(0) 

1 

r(°)  u  ) 

[U3  ’ 

Ul; 

[Uk  *  UlJ 

(1) 

> 

fd)  u  1 

b  * 

U2j 

k  2J 

u<2> 

3 

fu(2)  u  1 

uk  ’  3 

\  / 

(k-1) 


(5.51) 


and  the  matrix  U  by 


..(0)  _  v  A  -  i  k 

Uj  Yj  9  3  J*  9  •  *  •  9  K 


U1  1  ' 


u 


(0) 


(i)  (i-1)  ,  (i-1)  , 

uj  =  uj  -(uj  ’  ui}  ui 


u.  =  ^ 

J  J 


i  =  1,  . . . ,  k  -  1 


j  =  i  +  1,  . . . ,  k 


(5.52) 


The  subscript  notation  refers  to  the  corresponding  column  vectors 


in  Y  ,  ,  and  U  [see  equation  (5.43)]. 
old  new 

The  particular  solution  v  is  turned  into  the  orthogonal  complement 

of  the  new  base  set  U  (i.e.  the  vector, v  ,,,  is  converted  into  a  vector 

new  old 

v  such  that  v  is  orthogonal  to  all  the  components  of  U  ) 
new  new  new 


v 

new 


v  ,  ,  -  U  w 
old  new  new 


(5.53a) 
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where  wi  the  i-th  component  of  the  vector  w  is  defined  to  be 


W.  *  V  - ,  u.  . 
1  old  1 


(5.53b) 


Here  Vn l A  is  the  transpose  of  v  and  u.  is  the  i-th  column  of  U 

UJ*U  oxa  i  new 

It  should  be  noted  that  the  above  operation  simply  makes  use  of  the 
Gram-Schmidt  procedure. 

Several  tests  have  appeared  in  the  literature  to  determine  whether 
an  orthonormalization  is  necessary  (Scott  and  Watts  [62  ]  Gersting  and 
Jankowski  [22]).  Scott  and  Watts  have  suggested  a  check  which  is  described 
below.  The  algorithm  has  been  implemented  in  the  code  developed  for 
this  thesis.  A  slightly  different  version  is  used  in  the  Suport  code 
[62  ]•  The  basic  idea  is  to  avoid  using  solutions  that  violate  a 
condition  which  implies  that  linear  dependence  may  not  be  present. 

The  conditions  chosen  are, 


J-li 


<  e 


i  -  i,  2,  .. 


(5.54) 


and 


new 


<  e 


old 


where  e  =  10  and  L  is  about  three-fourths  of  the  number  of  significant 
decimal  digits  available  in  the  computer  wordlength.  Scott  and  Watts 
[62  ]  present  the  following  justification  of  the  above  criteria.  Suppose 
that  a  vector  in  the  old  set  is  dependent,  then  the  orthonormalization 

will  turn  it  into  a  null  vector •  This  provides  the  motivation 
for  comparing  the  size  of  P_  with  that  of  P^,  where  P^  =  ||  u°  ||  is 


(5.55) 
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used  since  the  vector  has  the  largest  magnitude  of  the  vectors  in 
the  old  set  as  a  result  of  the  pivoting  process. 

Since  the  magnitude  of  the  vectors  in  the  original  set  may  vary 
these  vectors  are  appropriately  scaled  by  producing  a  set  of  unit  vectors, 
YS.  Hence 

YS  =  UPS  =  UP  (5.56) 

where  S  =  diag  {  1/  ||  ||  ,  . . . ,  1/  ||  ||  } 


and 


P  =  Pjj—  .  (5.57) 


The  linear  dependence  test  on  the  diagonal  now  requires  that 


<  e  — ii-  =  e  .  (5.58) 

II  yj  II  II  ^ill 

The  above  expression  is  the  required  linear  dependence  test. 

(d)  Numerical  Solution  of  the  Mean  Energy  and  Mean  Dissipation  Equations 

The  system  of  equations  (3.67)— (3.71)  describing  the  mean  energy 
and  mean  dissipation  can  be  written  in  the  form 

+  G-  (K,  7,  y)  -  0  (5.59) 

is. 

+  G-  (K,  e,  y)  =0 


— 

S' 

dK 

• 

1  +  - 

ve 

dy 

J- 

r  -  ■« 

V 

1  +  t 

de 

va 

£ 

V  J 

dy 

_d_ 

dy 


(5.60) 
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These  differential  equations  are  subject  to  the  boundary  conditions 

K  =  0 ,  e  =  0  at  y  =  0  (5.61a) 

K  =  “t;  ,  e  =  “ r  for  large  y  (5.61b) 

C  ^  y 


where  k  is  the  Von  Karman  constant.  The  boundary  conditions  at  large  y 
are  simply  that  the  mean  flow  is  described  by  a  constant  stress  layer. 

In  the  above  equations  Gtt  and  G—  appear  as  source  terms  and  represent 
the  combined  effects  of  production  and  dissipation. 

These  highly  nonlinear  equations  are  solved  using  finite  difference 
techniques  and  Newton's  method.  In  order  to  derive  the  discretization 
equations  a  control  volume  is  set  up  on  a  variable  step  grid  such  that 
the  boundaries  are  midway  between  the  grid  points  (see  Figure  5.1). 
Equations  (5.59)  and  (5.60)  can  then  be  integrated  across  the  control 
volume  to  yield 
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The  derivatives  in  equations  (5.62)  and  (5.63)  are  approximated  by  using 

a  piecewise  linear  profile.  The  integral  terms  are  estimated  by  using 

the  average  value  of  the  functions  G  and  G_  respectively. 
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Typical  Control  Volume 


Figure  5.1.  Grid  Used  for  Integrating  Mean  K-e 
Model  Equations 
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Since  the  grid  consists  of  m  +  2  levels  and  since  the  energy  and 
dissipation  are  known  at  both  boundaries,  finite  difference  equations 
need  only  be  written  for  the  grid  points  2  ^  i  ^  m  +  1  , 
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The  above  system  represents  2  m  nonlinear  algebraic  equations  in 
2  m  unknowns  K  2  ,  ,  K3  ,  Km+2>  em+2’  The  al8ebraic  equations 

were  solved  using  NewtonTs  method,  Abrams  [  1  ].  The  solution  procedure 
was  to  linearize  the 'Source  terms  G  and  G  as  follows 
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and  to  linearize  the  diffusion  terms  as  follows, 
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Similar  expressions  to  (5.65a)  and  (5.65b)  can  be  derived  for  the 
diffusion  terms  (1  +  v^v)^  (K±n  -  K^)  and  (1  +  ^e±  ~ 

In  the  above  equations  the  superscript,  n,  refers  to  the  nth  iterate 
while  the  superscript,  n-1,  refers  to  the  previously  known  iterate.  For 
convenience  the  following  notation  is  introduced 


2  <  i  <  m  +  1  , 


(5.66a) 
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6_  =  S_  =  0  ,  (5.66c) 

K1  E1 

6_  =  6_  =  0.  (5 . 66d) 

^m+2  £m+2 

If  equations  (5 . 66a)- (5 . 66d)  are  substituted  into  equations 
(5 . 64a)- (5 . 65b)  and  the  resulting  equations  substituted  into  equations 
(5, 63a)- (5. 63d)  a  linear  system  of  2 m  equations  for  the  unknown  iterates 

d—  ,  6—  ,  . . .  ,  6—  .6—  *  is  obtained.  These  iterates  are 

2  2  m+1  nrt*l  — n— 1  — n— 1  — n— 1  — n— 1 

then  used  to  update  the  values  of  K2  ,  ^  *  •••»  ,  Gm-fl  ‘ 

This  procedure  is  repeated  until  a  satisfactory  level  of  convergence 


is  obtained. 


114 


CHAPTER  6 
RESULTS 

In  this  chapter  the  results  of  the  calculations  made  with  the 
various  models  developed  in  Chapter  3  are  presented.  An  interpretation 
and  discussion  of  the  results  is  made  in  the  final  chapter  of  this 
thesis . 

The  solution  of  equation  (3.13)  provides  detailed  information  on 
the  velocity,  shear  stress  and  pressure  fields.  It  is  found  that  the 
flow  field  in  the  immediate  vicinity  of  the  wave  surface  and  the  wave 
induced  surface  shear  stress  are  particularly  sensitive  to  the  type  of 
turbulence  model  employed.  Therefore,  emphasis  will  be  placed  on  the 
prediction  of  the  measured  surface  shear  stress  profiles  presented  in 
Chapter  4.  These  profiles  are  characterized  by  a  dimensionless  amplitude 
a  | x (o ) |  and  phase  angle  0  such  that  Real  (ax  (o)  elax)  =  a|x(o)|  cos 
(ax  +  0).  For  thick  boundary  layers  both  |x(o)|  and  0  are  functions  of 
a  single  parameter,  a,  the  dimensionless  wave  number. 

Where  possible  the  predicted  wave  induced  surface  pressure  is  also 
compared  with  the  available  measurements.  However,  the  considerable 
scatter  in  the  pressure  measurements  makes  the  comparison  difficult.  The 
wave  induced  pressure  is  also  characterized  by  a  dimensionless  amplitude 
a  | p (o ) |  and  a  phase  angle  0^  such  that  Real  (a p (o)  e^aX)  =  aJp(o)|  cos 
(ax  +  0  ). 

(a)  Quasi-Laminar  Model:  Model  A 

The  simplest  model  developed  in  Chapter  3  is  the  Quasi-Laminar 
Model  (Model  A),  for  which  the  effects  of  turbulence  enter  the  problem 
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only  in  the  specification  of  the  mean  velocity  profile.  Thorsness  [70] 
has  examined  several  velocity  profile  relations  and  concluded  that  the 
calculations  are  insensitive  to  the  choice  of  a  velocity  profile  function 
provided  it  gives  a  reasonable  fit  to  the  experimental  measurements. 
Consequently  the  calculations  were  carried  out  using  the  Van  Driest  velocity 
profile  with  the  Von  Karman  constant,  k  =  0.41,  and  the  viscous  thickness 
parameter,  A  =  26.  Figures  6.1  and  6.2  compare  the  predicted  values  of  the 
magnitude  and  phase  angle  of  the  wave  induced  shear  stress  with  measurements 
as  a  function  of  ct.  For  large  values  of  ct,  corresponding  to  small  wave¬ 
lengths,  experiment  and  theory  are  in  fair  agreement;  however, as  the  wave 
number  decreases  the  deviations  between  experiment  and  theory  increase. 

This  is  expected  since  the  influence  of  the  wave  surface  extends  further 
into  the  flow  as  the  wavelength  increases  (decreasing  a) .  For  large  values 
of  a  the  effects  of  the  wave  surface  are  confined  to  the  viscous  wall 
region.  Since  the  fully  turbulent  part  of  the  boundary  layer  is  undisturbed 
the  effect  of  turbulence  on  the  wave  induced  flow  is  not  important. 

Figures  6.3,  6.4,  6.5  and  6.6  show  the  wave  induced  velocity  profile 
function.  Real  (u(y)  eiaX) ,  for  a  =  0.1,  0.01,  0.0045  and  0.001  respectively. 
These  show  how  the  disturbance  caused  by  the  wave  is  confined  closer  to 
the  surface  as  a  increases. 

Figures  6.7  and  6.8  compare  the  predicted  magnitude  and  phase  angle 
of  the  surface  pressure  with  the  available  data. 

i 

(b)  Model  C* 

The  first  turbulence  model  tested.  Model  C*,  employs  the  Van  Driest 
mixing  length  hypothesis  to  estimate  a  turbulent  viscosity  and  ultimately 
a  turbulent  shear  stress.  The  model  uses  the  flat  plate  value  of  A  =  26 
for  the  viscous  sublayer  parameter.  To  account  for  the  effects  of  pressure 
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Figure  6.1.  |  t  (o) |  for  Quasi-Laminar  Model 
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Figure  6,7.  p (o)  |  for  Quasi-Laminar  Model 
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Figure  6.8.  p(o)  Phase  Angle  for  Quasi-Laminar  Model 
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gradient  on  the  thickness  of  the  viscosity  dominated  sublayer  the 
local  value  of  the  wall  shear  stress  is  used  in  the  Van  Driest  damping 
function.  Figures  6.9  and  6.10  compare  the  amplitude  and  phase  angle 
of  the  shear  stress  with  measured  values.  The  model  predicts  a  maximum 
in  the  phase  angle  of  the  shear  stress  at  a  =  0.007  which  is  not  what 
is  observed  experimentally.  The  model  also  fails  to  predict  the  sharp 
decrease  in  the  phase  angle  of  the  wave  induced  shear  stress. 

The  phase  angle  and  amplitude  of  the  wave  induced  surface  pressure 
are  shown  in  Figures  6.11  and  6.12.  Figures  6.13,  6.14  and  6.15  show 
the  wave  induced  velocity  profiles  for  a  =  0.01,  0.0045  and  0.001 
respectively. 

(c)  Equilibrium  Turbulence  Model 

The  next  model  tested  is  similar  to  Model  C*  except  that  the  effect 
of  pressure  gradient  on  the  viscous  sublayer  is  accounted  for  by  allowing 
the  sublayer  parameter.  A,  to  be  composed  of  a  mean  value  and  a  wave 
induced  valve,  AeiaX,  which  is  a  function  of  the  pressure  gradient.  This 
formulation  is  referred  to  as  the  Equilibrium  Turbulence  Model  since  the 
model  assumes  that  the  turbulence  adjusts  instantaneously  to  changes  in 
the  pressure  gradient. 

The  implementation  of  this  model  requires  the  specification  of  a 
single  constant,  defined  in  equation  3.47.  Kays  et  al .  suggest  a 

value  of  -30  <  <  “20  be  use<3.  Figure  6.16  shows  the  effect  of  this 

constant  on  the  calculated  phase  angle  of  the  wave  induced  shear  stress. 

For  an  equilibrium  situation  the  pressure  gradient  would  tend  to 
Shift  the  maximum  in  the  wall  shear  stress  in  the  downstream  direction 
because  an  adverse  pressure  gradient  which  tends  to  enhance  turbulence 
exists  on  the  downstream  side  of  the  wave.  Therefore  for  small  values 
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Figure  6,10,  t(o)  Phase  Angle  for  Model  C 
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of  Ot  the  consequence  of  including  pressure  gradient  effects  in  an 
equilibrium  model  is  to  give  lower  values  of  the  phase  shift  than  would 
be  obtained  with  Model  C*(k^  =  0),  Figure  6.16. 

(d)  Model  D*  and  Relaxation  Turbulence  Model 

The  turbulence  is  not  expected  to  be  in  a  state  of  equilibrium  with 
the  local  pressure  gradient.  Two  approaches  which  are  extensions  of  the 
Equilibrium  Turbulence  Model  are  explored. 

The  first  assumes  that  the  turbulence  only  in  the  viscous  wall  region 
is  not  in  equilibrium  with  the  local  pressure  gradient.  To  account  for 
this  an  effective  wave  induced  sublayer  parameter,  &  A  f  e  ax  ,  defined 
in  equation  (3.55)  is  used.  This  formulation  has  been  designated  as 
Model  D*  (see  Chapter  3). 

The  implementation  of  this  model  requires  the  specification  of  an 

additional  constant,  k^,  defined  in  equation  (3.53).  Figures  6.17  and 

6.18  compare  the  predicted  phase  angle  and  amplitude  of  the  wave  induced 

shear  stress  with  measurements.  Figures  6.19  and  6.20  compare  the 

predicted  phase  angle  and  amplitude  of  the  pressure  with  measurements. 

The  constants  used  in  these  calculations  are  k^  =  -35  and  k^  =  1800. 

Figures  6.21-6.24  show  the  effect  of  a  on  the  wave  induced  velocity  field. 

The  values  of  a  are  0.01,  0.0045,  0.001  and  0.0006  respectively. 

The  second  approach  differs  from  Model  D*  in  that  the  turbulence 

throughout  the  boundary  layer  is  not  considered  to  be  in  equilibrium  with 

the  pressure  gradient.  To  account  for  this  the  complete  Reynolds  stress 

a  iotx 

term  defined  in  equation  (3.33  )  is  relaxed,  rather  than  only  Ae 
The  implementation  of  this  model,  previously  referred  to  as  the  Relaxation 
Model,  also  requires  the  specification  of  a  second  constant,  kR,  defined 
in  equation  (3.14). 
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Figure  6.19.  p(o)  Phase  Angle  for  Model  D*  for  k 
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Figure  6.20.  |  p  (o) |  for  Model  D*  for  k-  =  -35  and  k  =  1800 
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An  excellent  fit  of  the  measurements  of  the  phase  angle  and  amplitude 

of  the  wave  induced  shear  stress  Figures  6.25  and  6.26  is  obtained  with 

k  =  -20  and  k„  =  4200.  Figures  6.27  and  6.28  show  the  effect  of  the 
1  R 

constants  and  ^  on  the  model  by  comparing  the  predicted  values  of  the 
fluctuating  wall  shear  stress  with  measurements  for  three  sets  of  values 
of  the  constants.  Figures  6.29  and  6.30  compare  the  predicted  phase 
angle  and  amplitude  of  the  wave  induced  pressure  with  measurements.  The 
effect  of  a  on  the  wave  induced  velocity  field  is  shown  in  Figures  6.31- 
6.34  for  a  »  0.01,  0.0045,  0.001  and  0.0006  respectively. 

(e)  K-e  Model 

The  last  model  to  be  tested  is  the  K-e  Model.  In  this  model  the 
effects  of  pressure  gradient  and  relaxation  are  accounted  for  in  a  natural 
way  through  the  defining  differential  equations.  The  mean  velocity  profile 
used  in  calculating  the  wave  induced  flow  is  not  the  Van  Driest  profile, 
but  rather  one  generated  using  the  mean  momentum,  mean  energy  and  mean 
dissipation  equations.  The  profiles  of  these  quantities  are  shown  in 
Figures  6.35-6.37  respectively. 

Figures  6.38  and  6.39  compare  the  calculated  phase  angle  and  amplitude 
of  the  wave  induced  shear  stress  predicted  by  the  K-e  Model.  The  calculated 
maximum  in  the  phase  angle  occurs  at  a  =  8  x  10  The  sharp  decrease  in 
the  phase  angle  as  a  function  of  a  is  not  as  apparent  as  that  predicted 
by  either  Model  D*  or  the  Relaxation  Model.  The  K-e  Model  also  under¬ 
predicts  the  magnitude  of  the  maximum  phase  shift  of  the  shear  stress. 
Figures  6.40  and  6.41  compare  the  predicted  phase  angle  and  amplitude  of 
the  wave  induced  pressure  with  measurements. 
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Figure  6.25.  x(o)  Phase  Angle  for  Relaxation  Model  for  k  =  -20 
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Figure  6.27,  Effect  of .  k  and  k  on  t(o)  Phase  Angle  for  Relaxation  Model 
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Figure  6.29,  p(o)  Phase  Angle  for  Relaxation  Model 
for  k.  =  -20  and  k  =  4200 
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Figure  6.30.  |  p(o)|  for  Relaxation  Model 
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Figure  6.32. 
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Figure  6,33. 
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Figure  6.34. 
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Figure  6.35.  Mean  Velocity  Profile  for  K-e  Model 
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Figure  6.36.  Mean  Energy  Profile  for  K-e  Model 
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Figure  6.37.  Mean  Dissipation  Profile  for  K-e  Model 
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Figure  6.38.  t(o)  Phase  Angle  for  K-e  Model 
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Figure  6*39.  x(o)|  for  K-e  Model 
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p(o)  Phase  Angle  for  K-e  Model 
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Figure  6.41.  |  p(o)  |  for  K-e  Model 
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In  Appendix  L  it  is  shown  how  the  wave  induced  turbulent  kinetic 
energy,  ak(y)elaX,  is  related  to  the  wave  induced  variation  of  the 
root  mean  square  level  of  the  fluctuating  velocity  gradient,  i/s^ /  Sx  . 
Figures  6.42  and  6.43  compare  the  phase  angle  and  amplitude  of  the 
predicted  wave  induced  kinetic  energy  with  the  experimentally  determined 
values.  In  the  range  of  the  experiments  the  predicted  amplitude  is  in 
fair  agreement  with  the  data.  Poor  agreement  between  the  measured  and 
predicted  phase  shifts  of  a  k  (y)  e^aX  is  obtained.  However  both  theory 
and  experiment  predict  a  phase  shift  such  that  the  maximum  occurs  on 
the  upstream  side  of  the  wave. 


(f )  Curvature  Effects 

A 

Model  C  was  modified  to  include  curvature  effects  by  multiplying 
the  plane  flow  mixing  length  expression  equation  (3.35)  by  a  curvature 
correction  equation  (3.79).  This  model  was  chosen  as  opposed  to  either 
the  Relaxation  Model  or  Model  D*  in  order  to  isolate  the  effect  of 
curvature  on  the  turbulence.  The  straightforward  application  of  the 
curvature  correction  requires  the  specification  of  one  constant,  3^, 
and  assumes  that  the  local  turbulence  is  always  in  equilibrium  with  the 
streamline  curvature. 

Figures  6.44,  6.45,  6.46  and  6.47  show  the  phase  Angle  and  amplitude 
of  the  i* ave  induced  pressure  and  shear  stress  for  different  values  of 
3^  respectively.  For  small  values  of  a  the  effect  of  curvature  Is  to 
shift  the  phase  angle  of  the  shear  stress  in  the  upstream  direction. 

This  may  be  explained  since  the  effect  of  the  wave  on  the  turbulence 
(from  curvature  considerations)  is  to  enhance  the  turbulence  at  the 
trough  and  suppress  the  turbulence  at  the  crest.  For  large  values  of 
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Figure  6.42.  k(o)  Phase  Angle  for  K-e  Model 
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Figure  6.43.  I  k(o)  |  for  K-e  Model 
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Figure  6.44.  p(o)  Phase  Angle  for  Model  C*  for  Different  Values  of 


Mode!  C  with  Curvature  Correction 
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Figure  6,46.  t(o)  Phase  Angle  for  Model  C*  for  Different  Values  of 
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Figure  6.47.  |  x(o)|  for  Model  C*  for  Different  Values  of 
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a  the  effect  of  curvature  is  to  shift  the  maximum  in  the  pressure  and 
hence  the  minimum  in  the  pressure  gradient  in  the  downstream  direction. 
This  in  turn  (from  laminar izat ion  arguments)  causes  the  shear  stress  to 
move  in  the  downstream  direction. 

Since  the  turbulence  is  not  expected  to  respond  instantaneously  to 
changes  in  the  streamline  curvature  Bradshaw  suggests  that  an  effective 
radius  of  curvature  defined  in  equation  (3.82)  be  used.  This  requires 
the  specification  of  a  curvature  lag  constant,^.  Figures  6.48,  6.49, 

6.50  and  6.51  show  the  effect  of  this  constant  on  the  phase  angle  and 
amplitude  of  the  wave  induced  pressure  and  shear  stress  respectively. 

Recent  pressure  measurements  from  Langley  Field  [35  ]  not  included 

in  this  thesis  indicate  that  the  phase  angle  measurements  of  Kendall 

are  the  most  reliable.  It  is  concluded  (see  Figure  6.48) that  inclusion 

of  curvature  with  (3^  *  1.5  provides  the  best  fit  of  the  available 

pressure  data.  Consequently  calculations  were  made  using  Model  D*  with 

a  curvature  correction.  This  model  was  chosen  since  it  allows  the 

turbulence  to  adjust  at  different  rates  to  the  effect  of  pressure  gradient 

and  streamline  curvature.  Figures  6.52,  6.53,  6.54  and  6.55  compare  the 

wave  induced  phase  angle  and  amplitude  of  the  shear  stress  and  pressure 

with  measurements.  The  constants  used  in  these  calculations  are  k.  32  -30, 

.I 

kT  a  1550,  3  =  2  and  k  * 

L  c  c 


0. 
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Figure  6.48.  p(o)  Phase  Angle  for  Model  C*  for  6=6  and  Different  Values  of  k 


Model  C*  with  Curvature  Correction 
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Figure  6,49.  |  p(o)|  for  Model  C*  for  3  =6  and 

Different  Values  k  C 


Model  C 
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Figure  6.50  t(o)  Phase  Angle  for  Model  C  for 
and  Different  Values  of  k 
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Figure  6.51.  |  t(o)|  for  Model  C*  for  6  =6  and  Different  Values  of  k 
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Figure  6.53*  |  t (o) |  for  Model  D*  for  k  -  -30,  =  1550, 
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Figure  6.54.  p(o)  Phase  Angle  for  Model  D*  for  k  =  -30,  k  =  1550, 


175 


M  rO 

o  c 

4~)  CO 

3  » 

<a  0 

-  GQ. 


176 


CHAPTER  7 

DISCUSSION  OF  RESULTS 

In  this  chapter  the  major  points  which  have  emerged  from  the 
experiments  and  the  calculations  are  presented. 

The  measurements  of  wall  shear  stress  have  been  shown  to  provide 
an  extremely  sensitive  test  for  turbulence  models.  Consequently  the 
reliability  of  any  theory  should  be  tested  by  examining  its  ability  to 
predict  the  variation  of  the  wall  shear  stress  over  a  wide  range  of 
conditions  rather  than  its  ability  to  predict  the  shear  stress  profiles 
for  some  fixed  flow  condition. 

Figures  7,1  and  7.2  compare  the  calculated  values  of  the  amplitude 
and  phase  angle  of  the  wave  induced  shear  stress  of  the  Quasi-Laminar 
Model  with  those  of  several  turbulence  models.  For  a  >  0.07,  the  Quasi- 
Laminar  Model  and  all  turbulence  models  yield  similar  results  since  the 

disturbance  field  introduced  by  the  wave  surface  is  contained  in  a  region 

•  #» 

where  the  effects  of  turbulence  on  the  wave  induced  flow  are  small. 

For  a  <  3  x  10  all  turbulence  models  predict  similar  variations  of  the 

shear  stress.  In  this  range  the  effect  of  pressure  gradient  on  the  flow 

is  small.  The  flow  can  be  visualized  as  one  over  a  surface  with  radius 

2 

of  curvature  1/R  »  -a  a  cos  (ax),  an  external  velocity  U  *  U  (I  +  a  cos  (ax)) 

2  2  s 

and  a  pressure  gradient,  3P/3x  =  -a  a  sin  (ax).  In  the  range  3  x  10 

<  a  <  0.03  the  Equilibrium  Turbulence  Model  over  predicts  the  effect  of 

pressure  gradient  and  the  only  models  which  reliably  predict  both  the 

magnitude  and  phase  angle  of  the  shear  stress  are  Model  D*  and  the 

Relaxation  Model  (not  shown  in  Figures  7.1  and  7.2).  Therefore  this  rise 

is  associated  with  relaxation  effects.  In  fact  for  a  >  7  x  10-^  relaxation 

effects  are  such  that  the  calculated  effective  pressure  gradient  can  be 
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Figure  7.1.  Comparison  of  Several  Calculation  Models  with  the  Data 
for  Shear  Stress  Phase  Angle 
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Figure  7.2.  Comparison  of  Several  Calculation  Models  with  the  Data  for  Shear  Stress  Amplitude 
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as  much  as  90°  out  of  phase  with  the  calculated  local  pressure  gradient 
(Figure  7.3)  causing  the  maximum  in  the  wall  shear  stress  to  shift  in 
the  upstream  direction. 

Measurements  of  the  streamwise  component  of  the  turbulent  kinetic 
energy  (see  Figure  6.42)  show  that  the  kinetic  energy  close  to  the  wave 
surface  reactes  a  maximum  on  the  upstream  side  of  the  wave.  This 
finding  provides  further  evidence  that  the  turbulence  close  to  the  wall 
in  the  range  3  x  10_5  <  a  <  0.02  is  not  in  equilibrium  with  the  local 
pressure  gradient,  since  (from  laminarization  arguments)  it  is  expected 
that  the  turbulent  energy  would  reach  a  maximum  on  the  downstream  side 
of  the  wave  due  to  flow  deceleration  in  this  region. 

The  inclusion  of  curvature  effects  in  Model  D  causes  a  significant 
improvement  in  the  calculated  phase  angle  of  the  wave  induced  pressure 
(Figure  7.4  in  the  range  1  x  10-4  <  a  <  0.5).  These  results  are  of 
particular  importance  in  the  prediction  of  instabilities  at  a  gas  liquid 
interface  where  recent  calculations  of  Frederick  [  21]  have  indicated 
that  Model  D*  may  be  underpredicting  the  pressure. 

In  summary  the  presence  of  waves  introduces  two  complications  that 

would  not  be  present  at  a  flat  surface:  Wave  induced  variations  of  the 

pressure  gradient  along  the  wave  surface  can  cause  a  thickening  and 

thinning  of  the  viscous  wall  region.  Wave  induced  curvature  of  the 

streamlines  can  cause  a  change  in  turbulent  transport  which  results  in 

a  periodic  enhancement  and  damping  of  the  turbulence.  Calculations 

have  shown  that  the  wall  shear  stress  measurements  and  in  particular 

,  -  .  .  iax 

the  sharp  change  in  the  measured  phase  angle  of,  ax  \o)  e  ,  at 
a  =  0.0015  provide  an  excellent  test  of  turbulence  models.  Significant 
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Figure  7.3.  Comparison  of  Phase  Angle  of  Pressure  Gradient  and  Effective 
Pressure  Gradient  for  Model  D* 
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errors  are  obtained  in  predicting  the  wall  shear  stress  if  the  influence 

of  the  wave  induced  pressure  gradient  is  ignored. 

The  sensitivity  of  Model  D*  to  the  value  of  the  relaxation  parameter, 

k  ,  points  out  the  need  to  introduce  relaxation  effects  in  a  more  funda- 
L 

mental  way.  Work  along  these  lines  using  the  K-e  Model  reveals  qualitatively 
the  same  behavior  as  Model  D*.  However  Model  D*,  modified  to  include  the 
effect  of  streamline  curvature  on  turbulence,  provides  the  best  overall 
fit  of  the  available  pressure  and  shear  stress  measurements  to  date. 
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APPENDIX  A 

CONDITIONS  FOR  LINEARITY 


The  necessary  conditions  for  flow  over  a  wavy  surface  to  be 
described  by  a  linear  form  of  the  Novier  Stokes  equations  are  developed. 

The  flow  field  is  considered  as  the  sum  of  a  mean  and  a  wave 
induced  quantity. 


U(x,  y) 

=  U(y)  +  u(x,  y) 

(A.  1) 

V(x,  y) 

=  v(x,  y) 

(A.  2) 

where  the  over  bar  and  tilda  represent  a  mean  and  a  wave  induced  quantity 
respectively.  The  wave  induced  streamwise  velocity  and  velocity  gradient 
scale  as 


dU  du 
dy  *  dx 


a  a 


dU 

dy 


(A.  3) 


since  the  presence  of  the  wave  surface  with  amplitude,  a,  displaces  the 
mean  flow  by  an  amount,  a  ^  ,  to  a  first  order  of  approximation.  From 
continuity 


9v  „ 
3y 


a  a 


dU 

dy 


(A.  4) 


On  integrating  (A.  4)  from  0  to  6,  where  6  is  a  boundary  layer  thickness 


V  ~ 


a  a  6 


dU 
dy  * 


(A.  5) 
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If  v  ~~  is  taken  as  a  typical  linear  term,  a  comparison  of  the 
ratio  of  nonlinear  to  linear  terms  yields 


3u 

3y 


-  dU 

v  dy 


a_ 

6 


(A.  6) 


~  3u 
U  3x 


dU 

dy 


_a 

6 


(A. 7) 


Therefore  in  addition  to  a  a  being  small,  the  additional  requirement  that 
a/6  <  1  is  necessary  for  the  flow  field  to  be  described  by  linear  equations. 
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APPENDIX  B 


FORMULATION  OF  OUTER  BOUNDARY  CONDITION 


The  outer  boundary  condition  is  such  that  the  flow  field  far  from 
the  surface  is  just  that  which  would  be  present  over  a  flat  surface. 
This  requires  that  the  velocity  field  in  boundary  layer  coordinates  at 
v  see  Figure  (B.l)be  equated  to  the  undisturbed  flow  in  cartesian 
coordinates. 

For  a  a  small  the  velocity  field  in  boundary  layer  coordinates  at 
ym  can  be  transformed  into  cartesian  coordinates  as  follows: 


VX’  ym}  =  U(x’  ym)  -  V(x’  ym} 


(B.l) 


(B.2) 


where  the  subscript,  o,  denotes  cartesian  coordinates.  If  U  and  V 
are  written  in  terms  of  average  and  wave  induced  quantities,  equations 


2 

(B.l)  and  (B.2)  neglecting  o(a  )  terms  become 


=  U(ym)  +  aGelotx 


(B.3) 


(B.4) 


The  cartesian  description  of  the  flow  at  (Y  +  aelax)  is  obtained  by  a 


Taylor  series  expansion  at  Y  as  follows 


m 


ae 


iax 


Y, 


m 


(B.5) 


V 


u 


V. 


Figure  B.l.  Diagram  for  Development  of  Outer  Boundary 
Condition 
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V 

o 


<x>  y«> 


=  V  (x,  Y  )  + 
o  m 


3Vq(x,  y) 


ae 


iax 


m 


(B.  6) 


Since  the  flow  at  Y  is  considered  undisturbed 

m 

U  (x,Y)  =  U(Y  )  (B. 7) 

o  in  m  v 

V  (x,  Y  )  =  0  -  (B.8) 

o  m 

Equations  (B.3)  and  (B.4)  can  now  be  equated  to  (B.5)  and  (B.6) 
respectively  to  yield 

U(y  )  +  au(y  )eictX  =  U(Y  )  +  ^  (Y  )ae±aX  (B.9) 

Jm  m  dy  m 

and 

U(y  )iaae^ax  +  vae*aX  -  0  .  (B,10) 

m 

On  equating  terms  of  similar  order  and  introducing  the  stream  function, 
iax 

aFe  ,  the  following  relations  are  obtained  , 


-  U(Y  ) 
m 

(B.ll) 

F'<y„) 

0 

>> 

|t3|  ^ 
T-d|tr) 

II 

(B.12) 

r<yn> 

/-\ 

it 

(B.13) 

Equation  (B.ll)  implies  that  the  average  flow  field  is  not  disturbed 
by  the  presence  of  the  wave  surface  and  equations  (B.12)  and  (B.13)  are 
the  required  boundary  conditions  for  the  wave  induced  flow  for  large 


values  of  y. 
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APPENDIX  C 


DERIVATION  OF  TURBULENT  ENERGY  AND  DISSIPATION  EQUATIONS 

Consider  a  flow  over  a  two  dimensional  curved  surface.  Use  general 
orthogonal  coordinates  with  x  measured  along  the  surface,  y  normal  to 
the  surface  and  z  at  right  angles  to  the  x-y  plane  which  is  the  plane 
of  motion.  The  curvature  of  the  surface  k(x)  =  1/R(x),  is  taken  as 
positive  for  convex  curvature  and  negative  for  concave  curvature  where, 
R(x) ,  is  the  radius  of  curvature  of  the  surface.  The  elements  of  length 

along  the  parallel  curves  and  along  the  normal  are  h^  =  1  +  <y  and 

h^  =  1.  The  element  of  length  along  the  z  direction  is  h^  ■  1.  If 
u,  v,  and  w  are  the  velocity  components  along  the  x,  y,  and  z  directions 

repectively  the  following  equations  for  u,  v,  and  w  are  obtained  (see 

Bradshaw  [  8  ] )  . 

continuity 


+  A  {(!  +  Ky)v}  +  A  {(i  +  xy)w}  =  0 


(c.l) 


x  momentum 


3u  1 

3t  (1  +  icy) 


3u 

3y 


(1  +  icy) 


K  UV 


1  3£ 

(1  +  icy)  3x 


(C.2) 
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y  momentum 


3v  1 

„  3v  4 

3t  (1  +  Ky) 

3x 

,  1 

3x 

xy 

(1  +  Ky) 

3x 

z  momentum 

3w  1 

3w 

3t  (1  +  <y) 

U  "r —  + 

3x 

i. 

3t 

xy  . 

(1  +  Ky) 

3x 

5v  3v  _  k  u 
3y  W  3z  (1  +  Ky) 


3y  3z  (1  +  <y)  xx  yy 


_  ®E 

3y 


(t _ -  x  _) 


3t  3t  3t  kt 


zz 


2L 


where 


x  =  2S  *  2 

XX  XX 


9u  +  <v 


(1  +  <y)  3x  (1  +  Ky) 


r  =  2S  =2  -|^ 
yy  yy  3y 


t  *  2S  -  2  ~ 
zz  zz  3z 


Txy  =  2Sxy  =  (1  + 


f  \ 

u 

l  1 

3v 

dy 

(1  +  Ky), 

(l  +  Ky) 

3x 

t  -  2S 
yz  yz 


3v  3w 
3z  3y  , 


t  =  2S 
xz  xz 


3u  1 


3w 


3z  (1  +  <y)  3xv 


(C.3) 


(C.4) 


(C.5) 


(C.6) 


(C.7) 


(C.8) 


(C.9) 


(C.IO) 


190 


Consider  a  turbulent  flow  that  is  two  dimensional  in  the  mean. 

Therefore,  the  mean  flow  in  the  z  direction  is  zero  and  all  —  of  the 

8z 

mean  flow  quantities  vanish.  If  the  equations  of  motion  are  decomposed 
into  mean  velocities  U,  V  plus  fluctuating  velocities  u1 ,  vT ,  w*  and 
their  corresponding  pressures  P  and  pf  the  following  equation  for  the 
turbulent  kinetic  energy  can  be  derived  (see  Bradshaw  [  8  ]). 

I 

9K  U  3K  ,  y  8K 

3t  (1  +  <y)  8x  3y 


II 


III 


1  8  f  uT  ^  +  uT vT  ^  -f  u TwT  ^  +  u*p1 

(1  +  <y)  8x  2 


Ufx !  -  V*X 1 

XX 


xy 


-  W T  T  1 
xz 


IV 


—  +  vT  ^  +  VTWf2  +  Vfpf 

2 


V 


UTI?  -VTTT  —  W 1 T 1  ) 

xy  yy  yz 


VI 


VII 


(1  +  Ky) 


^v'uf^  +  V*^  +  v T  w 1 ^  +  vfpf 


U1!1 


xy 


yy 


-  w'x' 


yz 


VIII 


81T 


K  U 


-  + 


3V 


.  9y  (1  +  Ky)  (1  +  Ky)  8x 


VIII 


VIII 


-u'2  T _ l _ + 

L  (1  +  Ky)  8x  (1  + 


kV _ 

Ky) 


-  v 


,2  JTV 
8y 
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IX 


xx 


3u'  ,  ,  3u' 

4  T  1  “ -  4 


3vf 


(1  4  icy)  3x  xy  3y  (1  4-  icy)  xy  3x 


IX 


,  f  3vf 

4-  T f  ~  4 


^ _  t  I 

yy  3y  (1  4-  <y)  xx 


(1  4-  icy)  “  xy 


IX 


TXZ  3w' 


(14-Ky)  3x  Tyz  3y  r  (1  +icy)  w  Lyz  Lzz  3z 


3w! 

f  -  4 


W  ?  T  *  -  t' 


3wf 


(C.ll) 


Jones  and  Launder  [29]  and  Chien  [13]  have  suggested  the  following 
assumptions 


3  2  2 

uf  4  ufvf  4  ufwf 


4  up 


t  3K 


(1  4  Ky)  va„  3x 
K. 


(C.12) 


t  f2  ,  t3  ,  f2 
rT  u  4  V  4  V  W 


4  v 1  p  1 


t  3K 


vaR  3y 


(C.13) 


The  terms  involving  the  fluctuating  shear  stress,  T±y  are  8enerally 
written  as  the  sum  of  a  term  representing  viscous  diffusion  of  kinetic 
energy  and  a  term  representing  the  isotropic  dissipation  of  kinetic  energy, 
Tennekes  and  Lumley  [69]. 


Ill  4  V  4  VII  4  IX  = 


4 


1 

1  3K  ^ 

+  J_ 

r  3k  ' 

(l  +  <y) 

3x 

[(1  +  <y)  3x  J 

3y 

(3y. 

K 

K  - 

2K 

(1  +  <y) 

e  2  * 
y 

(C.14) 
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2k 

where  the  terra  — ~  is  included  in  order  to  balance  the  viscous  diffusion 

y 

of  turbulent  kinetic  energy  y  at  the  wall. 

Equations  for  the  velocity  fluctuations  u* ,  v*  and  w1  can  be  written 
in  cartesian  tensor  notation  as  follows  (Tennekes  and  Lumley  [69]). 


3U 


i  +  ,  3Ui 

3t”  +  - 


au. 


3xk 


+  u; 


k  3x, 


+  u/ 


3u^  3u' 

3u! 

,  l 

3p' 

32ui 

*  i 

i 

c 

3xi 

3u! 

3u! 

c  dissipation 

rate  — ± 

3x. 

i 

3x . 

can  be 

(C. 15) 


obtained  by  differentiating  (C.15)  with  respect  to  x  ,  multiplying  the 

K 

resulting  equation  by  — —  and  taking  the  time  average 

3x 

n 


3c  ^  3e  ^  3 

3t  k  9x,  3xt 

K  K 


II 

3e 

3xk 


hi 


-  u;  e 


.jElK 

3x  3x 
n  n 


IV 


au. 


-  2 


)xk 


3u!  3u,' 

_ i  k 

3x  3x 
n  n 


3uf  3u! 
_ n  _ n 

3x±  3xk 


VI 


VII 


32U , 


3u! 


3u'  3 


-  2 


3x  “  3x.  3y  3x 
n  k  J  n  n 


-  2 


-  2 


•_3^ 

3VXn 


-T 


(C. 16) 


Jones  and  Launder  [29]  and  Chien  [13]  have  suggested  the  following 


closure  assumptions, 


193 


Term  III:  Turbulent  diffusion  by  velocity  and  pressure  fluctuations, 


III  = 


9  t  9e 


9xK  VCTe  9XK  ’ 


(C.17) 


Terms  V  and  VI:  Generation  by  mean  motion, 


V  +  VI  » 


ci£Pk 

K 


(C. 18) 


where  P  is  the  production  of  turbulent  kinetic  energy, 
K 


Term  VII:  Dissipation  of  e, 


e2  2e 

VII  =  C  f  - - —  f 

2  o  K  2  d  * 

y 


(C. 19) 


where  the  second  term  is  included  to  balance  viscous  diffusion  of  e  at 

the  wall.  The  functions  f  and  f ,  are  defined  in  Chapter  3  [see  equations 

o  d 

(3.63.1)  and  (3.63.2)]. 


Term  VI:  Generation  by  self  stretching  action  of  turbulence.  The  erfect 
of  this  term  is  accounted  for  in  the  modeling  of  term  IV  (see  Reynolds 
[59]). 

The  dissipation  equation  with  closure  assumptions  is  then  transformed 
into  boundary  layer  coordinates,  euqation  (3.62). 
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APPENDIX  D 

DERIVATION  OF  STREAMLINE  CURVATURE 


In  this  appendix  an  expression  for  the  curvature  of  the  streamlines 
in  boundary  layer  coordinates  is  derived.  Let  b(x)  be  a  position  vector 
describing  the  streamline  ij>(x,  y)  =  c  where 


<Kx,  y) 


•y  icXX 

U  (y)dy  +  F  ae  aX 
o 


(D.l) 


b(s)  =  ( x (s) ,  y(s)) 

If  s  is  a  measure  of  the  arc  length,  then  “  is  a  unit  vector  and  the 

magnitude  df  the  curvature  of  the  streamlines  i^(x,  y)  *  C  is  given  by 

d2S(s) 
ds2  ' 

In  boundary  layer  coordinates 


d| 

ds 


h 

x 


+ 


h  y'  1 

y  y 


(D.3) 


where  ?  and  ^  are  unit  vectors  along  and  perpendicular  to  the  wave  surface 
x  y 

respectively.  The  metrics  hx  and  h^  are  defined  as 


h  =  1  +  <y 
x 


(D.4.1) 


and 

h  *  1  o> 

y 

where  k  is  the  curvature  of  the  wave  surface.  The  primes  denote  differ¬ 
entiation  with  respect  to  s,  such  that 


,  dx(s) 

x  =  — - 

ds 


(D.5.1) 
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and 

v»  =  dyCsj 
y  ds 


(D.5.2) 


Therefore 


=  h1  xf  i"  +  h  x"  ix  +  xf  +  hf  yf  +  hy  y"  iy  +  h^  y1 
ds 


Also 


3 1 


31 


f<  m  x»  _2L  +  v 
x  x  3x  y  3y 


x 


31  31 

*  s  v*  ' '  —  -f-  v  *  - ^ 

y  X  3x  +  y  3y 


9ix  t  liB(1 

3x  K  1y  ’  3y 


81  t  3iy 

Ihf  “  K  x  ’  3y  "  ° 


Hence  (D.7.1)  and  (D.7.2)  become 


1'  =  -Kx'i 
x  y 


1'  =  KX’I 

y .  x 


Equation  (D.6)  can  now  be  written  as 


,2  * 


=  (h  x"  +  h'  x'  +  h  y'  k)  t  +  (h  y"  -  h  k(x')  ) 
,  z  x  x  y  x  y 

ds 


x-  ’>H 


(D.6) 

(D.7.1) 

(D.7.2) 

(D.8.1) 

(D.8.2) 

(D.9.1) 

(D.9.2) 


(D.10) 
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where  to  0(a) 

,  »  .3  iax  .  ,  2  iax  .  /TN  -  -  -  N 

=  aia  ye  x  +  aa  y  e  y  (D.ll.l) 

and 

h1  =  0  (D.11.2) 

y 

Since  ^  (x,  y)  =  C 

=  ijj  x*  +  \p  yf  =  0  .  (D.12) 

ds  x  y 

On  differentiating  (D.12)  with  respect  to  s  and  rearranging,  x"  and  y" 
can  be  written  explicitly  in  terms  of  x*  and  yT 


ip  x"  4*  ip  yn 
x  y 

38  - 

(x')2  -  ip  (y')2  -  2\j)  x' y'  , 
xx  yy  xy  J  ' 

(D.13.1) 

*  G. 

(D.13.2) 

db 

ds 

Since 

*  1  and  using  equation  (D.3) 

(hxx')2  +  (y*)2  =  1.  (D.  14) 


Now  on  differentiating  (D.14)  with  respect  to  s 

h^x!  x"  +  yf  y"  =  -h^x1  *  (D.15) 

(D.13.2)  and  (D.15)  can  now  be  solved  simultaneously  for  x"  and  y". 
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Gy'  -  *  (-h'  x') 
x"  =  - 2 - 5 - 

(D.16) 

ll;  y1  -  ip  h  x’ 

X  y  x 

ib  h'  x'  -  Gx' 
ii  x  X 

y  * 

(D.17) 

ip  y1  -  \p  h  y* 
x  y  x 

Now  to  0(a) 

\p  -  Fiaae^*aX 

X 

(D. 18) 

.  _  2  iax 

\p  =  -Fa  ae 

XX 

(D. 19) 

,  77  ,  iax 

ih  *  U  +  F  ae 

y 

(D. 20) 

.  tt  t  i  -nti  iax 

ip  =  UT  +  F"  ae 

yy 

(D. 21) 

#  „f  iax 

jp  ■  iaF1  ae 

Yxy 

(D. 22) 

.  ; ,  F  iax 

ip  /V  =  —  iaae 
x  yy  u 

(D. 22) 

Substituting  (D. 18)- (D. 23)  into  (D.12)  and  (D.14) 

and  neglecting 

2 

terms  of  0(a  )  xf  and  yT  can  be  expressed  in  terms  of 

aF  exax,  as  follows: 

the  stream  function, 

t  /i  2  i  (xx. 

x  *  (1  -  aa  ye  ) 

(D*  24) 

i  F  ,  iax 

y  *  iaa  e 

(D.25) 

U 


Substituting  (D. 18)- (D. 23)  into  (D.16)  and  (D.17)  and  neglecting  terms 
of  O(a^), 
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x"  =  iaya  e 


iax 


(D. 26) 


„  F  2  iax 
y  -  77  a  ae 


(D. 27) 


d^b 

Substituting  (D. 24)- (D. 27)  into  the  expression  for  — j  and  keeping  terms 

ds 

of  0(a)  the  following  equation  for  the  magnitude  of  the  curvature  is 
obtained 


d2b 


ds 


F  -  U 

2  iax 

aa  e 

U  J 

(D. 28) 


In  this  study  the  sign  of  the  curvature  is  defined  such  that  a 
convex  surface  has  a  positive  radius  of  curvature  and  a  concave  surface 
has  a  negative  radius  of  curvature.  With  this  convention 


_1_ 

R 

c 


(F  -  U)  2  iax 
- — — — -  aa  e 

U 


At  the  wave  surface  F  — ►  0  and  therefore 


(D. 29) 


R 

c 


2  iax 
aa  e 


(D. 30) 


which  is  simply  the  curvature  of  the  wave  surface.  For  large  y 

F  — *  U  and  here  — *  0 

K 

C 
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APPENDIX  E 

FINITE  BOUNDARY  LAYER  CALCULATION 


In  this  appendix  the  details  of  the  finite  boundary  layer  calculation 
are  presented.  As  explained  in  Chapter  3  a  finite  boundary  layer  consists  of 
two  layers.  In  the  inner  layer  the  mixing  length  is  given  by  equation 
(3.35)  and  in  the  outer  layer  the  mixing  length  is  given  by  equation  (3.90). 
In  order  to  provide  for  a  smooth  transition  between  these  two  layers  the 
following  composite  expression  is  used  across  the  entire  boundary  layer 


I/O  _ 

£  =  yd  tanh  (Ky/yd) (1  -  exp(-yx  /A)), 


(F.l) 


The  expression  for  £q  has  the  correct  asymptotic  behavior  since 


yd  tanh  (<y/y 6) 


<y 


y  < 


yd 


(E.2) 


and 


yd  tanh  (<y/yd)  — *■  yd 


yd 


(E.3) 


The  wave  induced  components  of  £  ,  v  and  t,  using  equations 
(E.l),  (3.34)  and  (3.36),  are  found  to  be  respectively 


a£  e 
o 


iax 


ayd  tanh  (<y/yd)  exp  (-y/A) 


_A 
A  j 


iax 


(E.  4) 


*  iax 
av^e  =  a 


s  2  exp  (-y/A) 

i  _ 

f  \ 

y 

'  ^ 

T 

A  ' 

,sxy  [1  -  exp (-y/A)] 

A 

w  2 

A  ^ 

_ 

v 

v 


—  eiaX  (E.5) 
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'  iax 
axe 


r  v 

2—  +  1 

V 


2  s 


xy 


2—  (uT  exp  (-y/A) 
[1  exp  (-y/A)  ] 


A 

A 


iax 


(E.6) 


1  - 


~  U'  exp  (— y / A)y 
[1  -  exp (-y/A)]  A 


For  large  values  of  y  equations  (E.4)  -  (E.6)  are  respectively 

a %  elax  =  0  (E.7) 

o 


(E.8) 

(E.9) 

Calculations  were  made  using  the  finite  boundary  layer  form  of 
Model  D  ,  equations  (E.4)-(E.6).  For  boundary  layer  thicknesses  6  <  500, 
it  is  not  possible  to  sustain  a  turbulent  boundary  layer  (see  Cebeci 
and  Smith  [12]).  Therefore,  this  was  the  smallest  value  of  6  considered. 

The  calculations  show  that  the  effect  of  6  on  the  wave  induced 
phase  angle  and  amplitude  of  the  shear  stress  and  pressure  is  small 
for  any  realistic  boundary  layer  thickness.  Figures  (E.l),  (E.2),  (E.3) 
and  (E.4)  respectively. 


t  iax  2^2.^  iax 

a  —  e  =  a  y  5  2  s  e 

v  xy 


-  iax  0  2.2^  /0  a  iax 

axe  =  a  2  y  <5  .U'  (2  s  )e 
1  xy 
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Figure  E.l.  Effect  of  Finite  Boundary  Layer  Thickness  on  t(o) 
Phase  Angle  for  Model  D*  k.  =  -35  and  =  1800 
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Figure  E.2.  Effect  of  Finite  Boundary  Layer  Thickness  on 
|  t  (o)  |  for  Model  D*  ^  =  -35  and  k  =  1800 
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(o)  d  jo  a|6uv  asoijd 


Figure  E.3.  Effect  of  Finite  Boundary  Layer  Thickness  on  p(o)  Phas 
Angle  for  Model  D*  k  =  -35  and  k  =  1800 
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Figure  E.4.  Effect  of  Finite  Boundary  Layer  Thickness  on 
I  P  (o ) I  for  Model  D*  k  =  -35  and  k  =  1800 

J-J 


APPENDIX  F 


TABULATED  DATA 

This  appendix  presents  in  tabulated  form  the  wave  profile, 
average  shear  stress  and  fluctuating  shear  stress  distributions 
reported  in  this  work. 

The  wave  profile  data  are  reported  in  inches.  The  shear 
stress  and  fluctuating  shear  stress  data  are  reported  in  the 
dimensionless  form  as  discussed  in  Chapter  4. 


MX 
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Table  F.l 


Wave  9  Wave  8 
in.  x  IQ3  in.  x  10 


0.00 

29.00 

30.00 

0.05 

28.00 

29.00 

0.10 

26.00 

28.00 

0.15 

24.00 

22.50 

0.20 

20.00 

22.00 

0.25 

16.00 

17.00 

0.30 

12.00 

13.00 

0.35 

8.00 

9.50 

0.40 

5.00 

6.00 

0.45 

1.50 

3.00 

0.50 

0.00 

0.00 

0.55 

0.50 

0.50 

0.60 

2.00 

3.00 

0.65 

6.00 

6.50 

0.70 

10.00 

10.00 

0.75 

14.00 

15.00 

0.80 

18.50 

19.00 

0.85 

23.00 

23.00 

0.90 

26.00 

27.00 

0.95 

28.00 

29.00 

1.00 

29.00 

30.00 

Wave  Profile  Measurements 

Wave  6 
in.  xlO3 

30.00 

28.50 
28.00 

22.50 
22.00 
17.00 
13.00 

9.  50 
6.00 
3.00 
0.00 
0.50 
3.00 

6.50 

10.50 
15.00 
19.00 

23.50 
27.00 
29.00 
30.00 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 
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Table  F.2.  Summary  of  Average  Shear  Stress  Measurements 


“b h/2 

2  It  v 

e 

T 

T  “  T 

W  W 

V 

,  * 

A  u  • 

~T 

W 

5970.00 

.00871 

56.82 

0.40 

5970.00 

.00871 

56.00 

0.40 

5970.00 

.00871 

56.70 

0.43 

6680.00 

.00787 

50.00 

0.38 

6680.00 

.00787 

54.50 

0.43 

6680.00 

.00787 

52.00 

0.39 

8450.00 

.00621 

57.20 

0.38 

8450.00 

.00621 

59.00 

0.42 

8450.00 

.00621 

60.00 

0.39 

8950.00 

.00607 

60.00 

0.36 

8950.00 

.00607 

63.40 

0.39 

8950.00 

.00607 

63.12 

0.36 

9650.00 

.00566 

60.30 

0.41 

9650.00 

.00566 

63.00 

0.41 

9650.00 

.00566 

56.50 

0.36 

11400.00 

.00489 

65.00 

0.32 

11400.00 

.00489 

65.50 

0.35 

11400.00 

.00489 

65.70 

0.34 

13000.00 

.00435 

55.46 

0.29 

13000.00 

.00435 

61.00 

0.34 

13000.00 

.00435 

64.14 

0.37 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 


Table  F.3.  Summary  of  Average  Shear  Stress  Measurements 


ub  h/2 

2  n  v 

V 

A  u 

14600.00 

.00389 

14600.00 

.00389 

14600.00 

.00389 

15700.00 

.00367 

15700.00 

.00367 

15700.00 

.00367 

19800.00 

.00297 

19800.00 

.00297 

19800.00 

.00297 

19800.00 

.00297 

19800.00 

.00297 

19800.00 

.00297 

20600.00 

.00287 

20600.00 

.00287 

20600.00 

.00287 

27950.00 

.00218 

27950.00 

.00218 

27950.00 

.00218 

30620.00 

.00201 

30620.00 

.00201 

30620.00 

.00201 

e 

T 

Tw"  Tw 

?w 

69.80 

0.31 

71.47 

0.34 

65.00 

0.36 

70.00 

0.31 

76.50 

0.33 

69.00 

0.32 

73.40 

0.27 

76.00 

0.27 

76.00 

0.25 

67.00 

0.29 

76.00 

0.29 

69.00 

0.24 

74.00 

0.23 

66.00 

0.29 

71.00 

0.27 

78.50 

0.20 

80.60 

0.21 

78.30 

0.20 

86.40 

0.16 

84.80 

0.18 

86.00 

0.17 

Table  F.4.  Summary  of  Average  Shear  Stress  Measurements 


Run  No. 

\  v 2 

2  nv 

9X 

Tw  Tw 

V 

,  * 

A  U 

Tw 

43 

34700.00 

.00181 

84.90 

0.19 

44 

34700.00 

.00181 

77.38 

0.20 

45 

36000.00 

.00174 

81.65 

0.16 

46 

36000.00 

.00174 

81.00 

0.14 

47 

36000.00 

.00174 

80.00 

0.16 

48 

36000.00 

.00174 

81.20 

0.17 

49 

36000.00 

.00174 

83.00 

0.16 

50 

58000.00 

.00135 

78.35 

0.13 

51 

58000.00 

.00135 

73.27 

0.15 

52 

58000.00 

.00135 

70.00 

0.12 

53 

58000.00 

.00135 

73.00 

0.11 

54 

64000.00 

.00139 

72.00 

0.23 

55 

64000.00 

.00139 

73.00 

0.11 

56 

71000.00 

.000951 

76.20 

0.11 

57 

71000.00 

.000951 

67.00 

0.13 

58 

71000.00 

.000951 

69.00 

0.09 

59 

73500.00 

.000912 

71.00 

0.10 

60 

73500.00 

.000912 

65.70 

0.08 

61 

73500.00 

.000912 

68.50 

0.09 

62 

73500.00 

.000912 

68.15 

0.08 

63 

84500.00 

.000813 

58.65 

0.08 

Table  F . 5.  Summary  of  Average  Shear  Stress  Measurements 


Run  No. 

nb  h/2 

2  n  v 

e 

T 

r  “  T 
w  w 

V 

.  * 

A  U 

aw 

64 

84500.00 

.000813 

54.10 

0.07 

65 

84500.00 

.000813 

57.00 

0.06 

66 

102500.00 

.000683 

45.00 

0.09 

67 

114500.00 

.000618 

28.00 

0.06 

68 

122500.00 

.000582 

17.00 

0.06 

69 

42500.00 

.00151 

79.00 

0.17 

70 

50000.00 

. .00129 

77.00 

0.13 

Table  F.6.  Summary  of  Fluctuating  Shear  Stress 
Measurements 


Run  No. 

ub  h/2 

2  n  v 

0 

/s'2  /  S 

X  X 

V 

Au* 

5 

6680.00 

.00787 

77.00 

.0523 

6 

6680.00 

.00787 

86.00 

.0402 

7 

8450.00 

.00621 

80.86 

.0448 

8 

8450.00 

.00621 

91.40 

.0497 

9 

8450.00 

.00621 

90.00 

.0501 

13 

9650.00 

.00566 

84.00 

.0502 

14 

9650.00 

.00566 

91.00 

.0560 

15 

9650.00 

.00566 

90.00 

.0351 

20 

13000.00 

.00435 

96.00 

.0316 

21 

13000.00 

.00435 

92.00 

.0416 

25 

15700.00 

.00367 

109.00 

.>0300 

26 

15700.00 

.00367 

119.00 

.0347 

27 

15700.00 

.00367 

108.00 

.0267 

32 

19800.00 

.00297 

117.00 

.0211 

35 

20600.00 

.00287 

100.00 

.0196 

36 

20600.00 

.00287 

119.00 

.0264 
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Table  F.7.  Average  Shear  Stress  Measurements 


X 

X 

Run  1 

Run  2 

Run  3 

Run  4 

Run  5 

Run  6 

0.00 

1.31 

1.20 

1.10 

1.17 

1.15 

1.21 

0.05 

1.12 

1.02 

1.21 

1.12 

1.11 

1.14 

0.10 

0.91 

0.99 

0.97 

1.05 

0.96 

1.00 

0.15 

0.82 

0.86 

0.87 

1.04 

0.92 

0.94 

0.20 

0.77 

0.75 

0.78 

0.76 

0.75 

0.78 

0.25 

0.66 

0.70 

0.68 

0.76 

0.71 

0.70 

0.30 

0.62 

0.62 

0.62 

0.61 

0.64 

0.62 

0.35 

0.61 

0.61 

0.61 

0.66 

0.60 

0.65 

0.40 

0.63 

0.63 

0.59 

0.59 

0.55 

0.60 

0.45 

0.65 

0.68 

0.63 

0.65 

0.61 

0.66 

0.50 

0.68 

0.77 

0.72 

0.68 

0.66 

0.72 

0.55 

0.85 

0.82 

0.79 

0.78 

0.78 

0.79 

0.60 

0.99 

0.93 

1.07 

0.94 

0.99 

0.88 

0.65 

1.15 

1.15 

1.12 

1.11 

1.18 

1.19 

0.70 

1.27 

1.28 

1.31 

1.20 

1.31 

1.24 

0.75 

1.35 

1.34 

1.32 

1.33 

1.42 

1.34 

0.80 

1.35 

1.40 

1.45 

1.34 

1.35 

1.39 

0.85 

1.40 

1.42 

1.32 

1.44 

1.42 

1.40 

0.90 

1.37 

1.34 

1.46 

1.35 

1.42 

1.29 

0.95 

1.31 

1.31 

1.23 

1.28 

1.32 

1.27 

1.00 

1.18 

1.22 

1.12 

1.14 

1.15 

1.21 
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Table  F.8.  Average  Shear  Stress  Measurements 


X 

X 

Run  7 

Run  8 

Run  9 

Run  10 

Run  11 

Run  12 

0.00 

1.26 

1.13 

1.10 

1.24 

1.14 

1.14 

0.05 

1.08 

1.08 

1.06 

1.10 

1.05 

1.07 

0.10 

0.95 

0.92 

0.91 

0.97 

0.92 

0.95 

0.15 

0.83 

0.85 

0.85 

0.85 

0.84 

0.85 

0.20 

0.70 

0.73 

0.76 

0.73 

0.74 

0.71 

0.25 

0.69 

0.65 

0.69 

0.70 

0.66 

0.66 

0.30 

0.61 

0.60 

0.69 

0.62 

0.63 

0.64 

0.35 

0.62 

0.64 

0.63 

0.67 

0.63 

0.66 

0.40 

0.65 

0.67 

0.66 

0.65 

0.62 

0.68 

0.45 

0.68 

0.70 

0.69 

0.70 

0.68 

0.72 

0.50 

0.76 

0.74 

0.72 

0.76 

0.77 

0.81 

0.55 

0.89 

0.81 

0.74 

0.91 

0.88 

0.88 

0.60 

1.03 

1.14 

0.98 

1.05 

1.10 

1.01 

0.65 

1.13 

1.09 

1.21 

1.20 

1.20 

1.20 

0.70 

1.25 

1.34 

1.34 

1.25 

1.32 

1.32 

0.75 

1.35 

1.41 

1.39 

1.33 

1.33 

1.34 

0.80 

1.35 

1.30 

1.46 

1.30 

1.39 

1.35 

0.85 

1.26 

1.41 

1.39 

1.33 

1.35 

1.33 

0.90 

1.38 

1.46 

1.34 

1.31 

1.36 

1.27 

0.95 

1.33 

1.15 

1.26 

1.21 

1.23 

1.24 

1.00 

1.21 

1.20 

1.15 

1.13 

1.15 

1.16 
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Table  F.9.  Average  Shear  Stress  Measurements 


X 

X 

Run  13 

Run  14 

Run  15 

Run  16 

Run  17 

Run  18 

0.00 

1.16 

1.12 

1.16 

1.20 

1.18 

1.09 

0.05 

1.06 

1.04 

1.08 

1.08 

1.02 

1.04 

0.10 

0.94 

0.91 

0.97 

0.95 

1.04 

0.93 

0.15 

0.81 

0.81 

0.88 

0.84 

0.83 

0.83 

0.20 

0.71 

0.74 

0.78 

0.72 

0.67 

0.69 

0.25 

0.67 

0.66 

0.68 

0.70 

0.67 

0.68 

0.30 

0.60 

0.64 

0.65 

0.66 

0.67 

0.67 

0.35 

0.63 

0.63 

0.65 

0.71 

0.66 

0.69 

0.40 

0.62 

0.62 

0.68 

0.71 

0.69 

0.72 

0.45 

0.66 

0.63 

0.72 

0.74 

0.71 

0.77 

0.50 

0.73 

0.76 

0.77 

0.83 

0.83 

0.83 

0.55 

0.88 

0.80 

0.81 

0.95 

0.97 

0.96 

0.60 

1.03 

1.07 

0.95 

1.09 

1.02 

1.04 

0.65 

1.17 

1.22 

1.14 

1.18 

1.20 

1.22 

0.70 

1.32 

1.37 

1.29 

1.27 

1.31 

1.30 

0.75 

1.37 

1.39 

1.30 

1.31 

1.33 

1.32 

0.80 

1.37 

1.45 

1.36 

1.25 

1.35 

1.31 

0.85 

'  1.39 

1.41 

1.34 

1.29 

1.29 

1.31 

0.90 

1.39 

1.34 

1.34 

1.26 

1.31 

1.28 

0.95 

1.30 

1.26 

1.27 

1.15 

1.18 

1.19 

1.00 

1.19 

1.14 

1.16 

1.11 

1.08 

1.13 
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Table  F.10.  Average  Shear  Stress  Measurements 


X 

X 

Run  19 

Run  20 

Run  21 

Run  22 

Run  23 

Run  24 

0.00 

1.22 

1.16 

1.14 

1.18 

1.07 

1.06 

0.05 

1.15 

1.09 

1.06 

1.02 

1.04 

1.02 

0.10 

1.00 

0.96 

0.94 

0.95 

0.90 

0.91 

0.15 

0.91 

0.85 

0.83 

0.79 

0.80 

0.81 

0.20 

0.78 

0.77 

0.75 

0.74 

0.76 

0.70 

0.25 

0.75 

0.69 

0.68 

0.68 

0.66 

0.68 

0.30 

0.71 

0.66 

0.67 

0.67 

0.68 

0.67 

0.35 

0.69 

0.69 

0.67 

0.70 

0.67 

0.72 

0.40 

0.71 

0.71 

0.65 

0.76 

0.69 

0.75 

0.45 

0.79 

0.74 

0.63 

0.81 

0.77 

0.78 

0.50 

0.81 

0.76 

0.75 

0.85 

0.87 

0.85 

0.55 

0.90 

0.85 

0.89 

0.96 

0.94 

0.95 

0.60 

1.03 

1.01 

1.12 

1.13 

1.15 

1.09 

0.65 

1.03 

1.20 

1.19 

1.16 

1.20 

1.23 

0.70 

1.26 

1.30 

1.35 

1.34 

1.33 

1.32 

0.75 

1.26 

1.31 

1.33 

1.25 

1.31 

1.30 

0.80 

1.20 

1.35 

1.38 

1.29 

1.34 

1.31 

0.85 

1.26 

1.30 

1.31 

1.24 

1.26 

1.30 

0.90 

1.22 

1.27 

1.37 

1.25 

1.29 

1.23 

0.95 

1.21 

1.20 

1.21 

1.15 

1.17 

1.18 

1.00 

1.10 

1.14 

1.07 

1.10 

1.10 

1.14 

>"  I  X 
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Tabl£  F.ll.  Average  Shear  Stress  Measurements 


t  /  T 

W  W 


Run  25  Run  26  Run  27  Run  28  Run  29  Run  30 


0.00 

0.05 

0.10 

0.15 

0.20 

0.25 

0.30 

0.35 

0.40 

0.45 

0.50 

0.55 

0.60 

0.65 

0.70 

0.75 

0.80 

0.85 

0.90 

0.95 

1.00 


1.10 

1.04 

0.92 

0.81 

0.73 

0.69 

0.67 

0.74 

0.74 

0.77 

0.83 

1.00 

1.11 

1.23 

1.25 

1.26 

1.26 

1.26 

1.26 
1.25 
1.07 


1.05 

0.98 

0.86 

0.81 

0.74 

0.69 

0.70 

0.72 

0.75 

0.80 

0.86 

0.98 

1.15 

1.23 

1.31 

1.34 

1.31 
1.28 

1.32 
1.18 
0.95 


1.09 

1.07 

0.94 

0.84 

0.73 

0.63 

0.68 

0.70 

0.78 

0.81 

0.84 

0.92 

1.09 

1.24 

1.32 

1.34 

1.23 

1.29 

1.16 

1.15 

1.14 


1.13 

1.02 

0.91 

0.81 

0.76 

0.71 

0.71 

0.75 

0.81 

0.84 

0.87 

1.01 

1.08 

1.19 

1.21 

1.28 

1.21 

1.27 

1.19 

1.18 

1.04 


1.08 

0.95 

0.93 

0.81 

0.77 

0.71 

0.74 

0.76 

0.81 

0.83 

0.93 

1.01 

1.12 

1.16 

1.25 

1.22 

1.29 

1.23 

1.21 

1.13 

1.07 


1.10 

1.00 

0.95 

0.86 

0.79 

0.71 

0.72 

0.74 

0.81 

0.87 

0.92 

0.97 

1.17 

1.15 

1.27 

1.30 

1.17 

1.26 

1.08 

1.06 

1.11 
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Table  F.12.  Average  Shear  Stress  Measurements 


X 

A 

Run  31 

Run  32 

Run  33 

Run  34 

Run  35 

Run  36 

0.00 

1.16 

1.05 

1.10 

1.12 

1.12 

1.14 

0.05 

1.02 

0.98 

1.04 

1.06 

1.07 

1.09 

0.10 

0.94 

0.88 

0.93 

0.94 

0.91 

0.93 

0.15 

0.80 

0.83 

0.88 

0.85 

0.88 

0.86 

0.20 

0.77 

0.76 

0.80 

0.76 

0.77 

0.78 

0.25 

0.70 

0.75 

0.74 

0.78 

0.74 

0.66 

0.30 

0.72 

0.75 

0.75 

0.75 

0.71 

0.69 

0.35 

0.74 

0.75 

0.81 

0.80 

0.76 

0.67 

0.40 

0.79 

0.75 

0.80 

0.83 

0.71 

0.86 

0.45 

0.80 

0.85 

0.83 

0.87 

0.75 

0.86 

0.50 

0.85 

0.90 

0.89 

0.87 

0.84 

0.88 

0.55 

0.95 

0.97 

0.95 

1.02 

0.94 

0.94 

0.60 

1.08 

1.13 

1.03 

1.08 

1.12 

1.10 

0.65 

1.14 

1.19 

1.22 

1.15 

1.17 

1.23 

0.70 

1.22 

1.25 

1.14 

1.20 

1.26 

1.25 

0.75 

1.25 

1.31 

1.24 

1.23 

1.28 

1.29 

0.80 

1.25 

1.28 

1.27 

1.18 

1.26 

1.21 

0.85 

1.25 

1.25 

1.25 

1.23 

1.19 

1.20 

0.90 

1.25 

1.29 

1.16 

1.18 

1.25 

1.10 

0.95 

1.21 

1.16 

1.09 

1.11 

1.14 

1.20 

1.00 

1.12 

0.93 

1.08 

0.98 

1.12 

1.07 
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Table  F.13.  Average  Shear  Stress  Measurements 

x  /t 
w  w 


Run  37  Run  38  Run  39  Run  40  Run  41  Run  L2 


x 

X 

0.00 

0.05 

0.10 

0.15 

0.20 

0.25 

0.30 

0.35 

0.40 

0.45 

0.50 

0.55 

0.60 

0.65 

0.70 

0.75 

0.80 

0.85 

0.90 

0.95 

1.00 


1.08 

1.00 

0.93 

0.84 

0.80 

0.81 

0.80 

0.86 

0.88 

0.89 

0.92 

1.00 

1.10 

1.14 

1.19 

1.19 

1.16 

1.16 

1.16 

1.08 

1.01 


1.03 

0.99 

0.89 

0.86 

0.81 

0.79 

0.79 

0.84 

0.87 

0.92 

0.95 

0.99 

1.11 

1.13 

1.21 

i.h 

1.19 

1.16 

1.18 

1.08 

1.00 


1.02 

0.98 

0.91 

0.85 

0.80 

0.79 

0.83 

0.87 

0.89 

0.84 

0.94 

0.97 

1.09 

1.16 

1.20 

1.18 

1.20 

1.17 

1.14 

1.10 

1.05 


1.0$ 

1.02 

0.90 

0.86 

0.81 

0.85 

0.80 

0.91 

0.88 

0.93 

0.-98 

1.04 

1.08 

1.14 
1.17 
1.19 
1.10 

1.15 
1.10 
1.03 
1.00 


0.99 

1.00 

0.88 

0.88 

0.82 

0.84 

0.81 

0.87 

0.90 

0.93 

0.96 

0.99 

1.10 

1.13 

1.18 

1.20 

1.18 

1.18 

1.16 

1.08 

0.92 


0.87 

1.03 

0.91 

0.8$ 

0.81 

0.84 

0.84 

0.92 

0.91 

0.92 

0.98 

1.02 

1.06 

1.12 

1.16 

1.18 

ia5 

1.17 

1.10 

1.08 

1.03 
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Table  F.14.  Average  Shear  Stress  Measurements 


X 

X 

Run  43 

Run  44 

Run  45 

Run  46 

Run  47 

Run  48 

0.00 

0.97 

1.01 

1.07 

1.00 

1.05 

1.02 

0.05 

0.98 

1.01 

1.01 

1.01 

0.97 

0.97 

0.10 

0.86 

0.89 

0.92 

0.91 

0.90 

0.91 

0.15 

0.89 

0.92 

0.87 

0.92 

0.85 

0.90 

0.20 

0.81 

0.87 

0.84 

0.85 

0.83 

0.86 

0.25 

0.87 

0.87 

0.85 

0.89 

0.84 

0.82 

0.30 

0.84 

0.85 

0.82 

0.86 

0.86 

0.86 

0.35 

0.86 

0.92 

0.90 

0.93 

0.90  ' 

0.85 

0.40 

0.87 

0.86 

0.91 

0.90 

0.90 

0.91 

0.45 

0.89 

0.88 

0.93 

0.93 

0.91 

0.96 

0.50 

0.91 

0.96 

0.96 

0.96 

0.95 

0.99 

0.55 

0.99 

0.98 

1.01 

1.00 

1.01 

0.97 

0.60 

1.07 

0.99 

1.08 

1.03 

1.06 

1.06 

0.65 

1.17 

1.18 

1.12 

1.14 

1.12 

1.07 

0.70 

1.14 

1.06 

1.15 

1.13 

1.15 

1.18 

0.75 

1.26 

1.19 

1.18 

1.16 

1.17 

1.19 

0.80 

1.16 

1.16 

1.14 

1.12 

1.14' 

1.14 

0.85 

1.21 

1.21 

1.13 

1.16 

■  1.15 

1.11 

0.90 

1.16 

1.08 

1.11 

1.08 

1.13 

1.19 

0.95 

1.12 

1.06 

1.06 

1.07 

1.09 

1.06 

1.00 

0.97 

1.05 

0.96 

1.00 

1.01 

0.99 
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Table  F.15.  Average  Shear  Stress  Measurements 


X 

A 

Run  49 

Run  50 

Run  51 

Run  52 

Run  53 

Run  54 

0.00 

1.01 

1.06 

1.01 

1.03 

1.02 

1.06 

0.05 

0.96 

1.00 

1.02 

1.03 

1.04 

1.00 

0.10 

0.91 

0.95 

0.92 

1.00 

0.98 

0.99 

0.15 

0.89 

0.87 

0.87 

0.92 

0.92 

0.92 

0.20 

0.86 

0.85 

0.90 

0.87 

0.90 

0.89 

0.25 

0.84 

0.92 

0.92 

0.89 

0.90 

0.95 

0.30 

0.86 

0.82 

0.91 

0.89 

0.91 

0.90 

0.35 

0.90 

0.92 

0.92 

0.93 

0.94 

0.91 

0.40 

0.90 

0.93 

,0.82 

0.91 

0.93 

0.94 

0.45 

0.93 

0.92 

0.84 

0.93 

0.93 

0.92 

0.50 

0.95 

0.95 

0.86 

0.94 

0.95 

0.96 

0.55 

0.98 

0.98 

0.91 

0.95 

0.96 

0.97 

0.60 

1.07 

1.04 

1.00 

0.96 

1.03 

1.03 

0.65 

1.07 

1.08 

1.07 

1.11 

1.08 

1.07 

0.70 

1.15 

1.16 

1.19 

1.12 

1.09 

1.12 

0.75 

1.15 

1.15 

1.18 

1.12 

1.12 

1.10 

0.80 

1.16 

1.08 

1.18 

1.12 

1.12 

1.04 

0.85 

1.16 

1.13 

1.19 

1.11 

1.11 

1.10 

0.90 

1.15 

1.11 

1.21 

1.07 

1.05 

1.10 

0.95 

1.09 

1.09 

1.12 

1.05 

1.01 

1.06 

1.00 

1.02 

0.98 

0.99 

1.05 

1.00 

0.98 
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Table  F.16.  Average  Shear  Stress  Measurements 


X 

X 

Run  55 

Run  56 

Run  57 

Run  58 

Run  59 

Run  60 

0.00 

0.99 

1.04 

1.02 

1.07 

1.04 

1.01 

0.05 

1.00 

0.98 

1  j04 

1.04 

1.01 

1.02 

0.10 

0.94 

0.93 

0.94 

1.02 

0.97 

0.97 

0.15 

0.95 

0.88 

0.89 

0.94 

0.92 

0.98 

0.20 

0.93 

0.90 

0.98 

0.89 

0.91 

0.97 

0.25 

0.92 

0.94 

0.93 

0.92 

0.90 

0.94 

0.30 

0.90 

0.85 

0.92 

0.90 

0.90 

0.92 

0.35 

0.95 

0.93 

0.84 

0.93 

0.94 

0.96 

0.40 

0.94 

0.95 

0.86 

0.89 

0.94 

0.93 

0.45 

0.92 

0.94 

0.87 

0.94 

0.94 

0.92 

0.50 

0.96 

0.91 

0.84 

0.96 

0.96 

0.93 

0.55 

0.97 

0.97 

0.87 

0.99 

0.96 

0.96 

0.60 

0.98 

1.01 

0.96 

1.05 

1.02 

0.99 

0.65 

1.09 

1.05 

1.05 

1.08 

1.04 

1.07 

0.70 

1.10 

1.16 

1.16 

1.08 

1.10 

1.08 

0.75 

1.11 

1.16 

1.16 

1.07 

1.08 

1.09 

0.80 

1.13 

1.07 

1.21 

1.07 

1.10 

1.11 

0.85 

1.09 

1.12 

1.16 

1.08 

1.10 

1.07 

0.90 

1.06 

1.11 

1.19 

1.03 

1.11 

1.05 

0.95 

1.04 

1.09 

1.09 

1.01 

1.07 

1.03 

1.00 

1.04 

1.00 

1.02 

1.02 

0.99 

1.02 
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Table  F.17.  Average  Shear  Stress  Measurements 


X 

X 

Run  61 

Run  62 

Run  63 

Run  64 

Run  65 

Run  66 

0.00 

1.06 

1.00 

1.03 

0.97 

0.98 

1.06 

0.05 

1.01 

1.01 

1.03 

1.01 

1.03 

1.01 

0.10 

0.99 

0.97 

0.98 

0.93 

0.99 

1.03 

0.15 

0.91 

0.96 

0.96 

0.93 

1.00 

0.97 

0.20 

0.92 

0.96 

0.91 

0.91 

0.97 

0.95 

0.25 

0.95 

0.95 

0.99 

0.97 

0.96 

0.94 

0.30 

0.90 

0.92 

0.89 

0.92 

0.94 

0.91 

0.35 

0.93 

0.96 

0.95 

0.94 

0.98 

0.93 

0.40 

0.95 

0.93 

0.91 

0.96 

0.93 

0.95 

0.45 

0.93 

0.92 

0.95 

0.91 

0.92 

0.91 

0.50 

0.95 

0.94 

0.96 

0.92 

0.93 

0.92 

0.55 

0.96 

0.96 

0.96 

0.97 

0.93 

0.94 

0.60 

1.02 

0.97 

0.97 

1.05 

0.93 

0.99 

0.65 

1.06 

1.07 

1.04 

1.09 

1.05 

1.01 

0.70 

1.10 

1.09 

1.07 

1.12 

1.04 

1.07 

0.75 

1.08 

1.09 

1.09 

1.12 

1.09 

1.04 

0.80 

1.03 

1.12 

1.00 

1.15 

1.09 

1.08 

0.85 

1.09 

1.07 

1.12 

1.15 

1.10 

1.10 

0.90 

1.10 

1.06 

1.09 

1.09 

1.04 

1.11 

0.95 

1.08 

1.05 

1.11 

0.97 

1.04 

1.10 

1.00 

0.99 

1.00 

0.98 

0.93 

1.06 

1.05 
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Table  F.18.  Average  Shear  Stress  Measurements 


X 

X 

Run  67 

Run  68 

0.00 

1.03 

1.04 

0.05 

1.08 

1.09 

0.10 

1.06 

1.06 

0.15 

1.02 

1.03 

0.20 

0.97 

0.97 

0.25 

1.01 

1.00 

0.30 

0.93 

0.93 

0.35 

0.96 

0.97 

0.40 

0.94 

0.94 

0.45 

0.96 

0.92 

0.50 

0.94 

0.94 

0.55 

0.93 

0.93 

0.60 

0.95 

0.96 

0.65 

0.98 

1.01 

0.70 

1.03 

1.03 

0.75 

1.05 

1.03 

0.80 

1.02 

1.03 

0.85 

1.03 

1.03 

0.90 

1.02 

1.02 

0.95 

1.04 

1.03 

1.00 

1.02 

1.01 
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Table  F.19.  Fluctuating  Shear  Stress  Measurements 


Run  5  Run  6  Run  7  Run  8  Run  9  Run  15 


x 

X 

0.00 

0.05 

0.10 

0.15 

0.20 

0.25 

0.30 

0.35 

0.40 

0.45 

0.50 

0.55 

0.60 

0.65 

0.70 

0.75 

0.80 

0.85 

0.90 

0.95 

1.00 


0.29 

0.29 

0.27 

0.30 

0.25 

0.27 

0.26 

0.25 

0.24 

0.27 

0.28 

0.28 

0.32 

0.36 

0.37 

0.37 

0.33 

0.35 

0.34 

0.33 

0.30 


0.31 

0.31 

0.29 

0.30 

0.27 

0.26 

0.26 

0.28 

0.26 

0.29 

0.31 

0.31 

0.32 

0.38 

0.35 

0.34 

0.34 

0.34 

0.31 

0.52 

0.31 


0.32 

0.30 

0.27 

0.27 

0.24 

0.26 

0.24 

0.26 

0.27 

0.29 

0.29 

0.31 

0.32 

0.33 

0.34 

0.36 

0.34 

0.30 

0.34 

0.34 

0.31 


0.28 

0.30 

0.26 

0.27 

0.24 

0.25 

0.24 

0.26 

0.28 

0.31 

0.28 

0.29 

0.34 

0.34 

0.36 

0.36 

0.31 

0.34 

0.36 

0.28 

0.31 


0.28 

0.29 

0.26 

0.26 

0.25 

0.25 

0.28 

0.27 

0.27 

0.29 

0.29 

0.28 

0.33 

0.36 

0.36 

0.36 

0.36 

0.34 

0.33 

0.32 

0.30 


0.30 

0.30 

0.29 

0.28 

0.27 

0.26 

0.26 

0.28 

0.28 

0.30 

0.31 

0.31 

0.31 

0.33 

0.34 

0.34 

0.34 

0.32 

0.32 

0.32 

0.29 
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Table  F.20.  Fluctuating  Shear  Stress  Measurements 
/s'2/S 

V  X  X 


X 

A 

Run  13 

Run  16 

Run  20 

Run  21 

Run  25 

Run  2 

0.00 

0.30 

0.29 

0.29 

0.28 

0.27 

0.25 

0.05 

0.29 

0.28 

0.29 

0.29 

0.29 

0.25 

0.10 

0.28 

0.27 

0.28 

0.27 

0.26 

0.25 

0.15 

0.26 

0.25 

0.27 

0.27 

0.25 

0.25 

0.20 

0.25 

0.26 

0.27 

0.26 

0.24 

0.24 

0.25 

0.25 

0.24 

0.26 

0.25 

0.25 

0.24 

0.30 

0.25 

0.25 

0.26 

0.26 

0.25 

0.25 

0.35 

0.27 

0.26 

0.27 

0.27 

0.28 

0.28 

0.40 

0.26 

0.27 

0.27 

0.26 

0.28 

0.28 

0.45 

0.28 

0.26 

0.28 

0.26 

0.29 

0.30 

0.50 

0.28 

0.31 

0.29 

0.28 

0.29 

0.28 

0.55 

0.31 

0.27 

0.31 

0.29 

0.30 

0.30 

0.60 

0.33 

0.34 

0.32 

0.34 

0.32 

0.32 

0.65 

0.35 

0.35 

0.33 

0.33 

0.31 

0.31 

0.70 

0.36 

0.38 

0.33 

0.35 

0.32 

0.31 

0.75 

0.35 

0.36 

0.32 

0.34 

0.31 

0.32 

0.80 

0.35 

0.35 

0.33 

0.33 

0.29 

0.30 

0.85 

0.35 

0.34 

0.31 

0.32 

0.29 

0.29 

0.90 

0.34 

0.33 

0.30 

0.33 

0.30 

0.30 

0.95 

0.33 

0.31 

0.30 

0.30 

0.30 

0.29 

1.00 

0.31 

0.29 

0.29 

0.27 

0.26 

0.24 
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table  F.21.  Fluctuating  Shear  Stress  Measurements 


X 

X 

Run  27 

Run  32 

Run  35 

Run  36 

0.00 

0.28 

0.25 

0.26 

0.26 

0.05 

0.28 

0.25 

0.26 

0.26 

0.10 

0.27 

0.23 

0.24 

0.25 

0.15 

0.25 

0.24 

0.25 

0.25 

0.20 

0.25 

0.24 

0.24 

0.25 

0.25 

0.23 

0.25 

0.24 

0.21 

0.30 

0.25 

0.26 

0.24 

0.24 

0.35 

0.26 

0.25 

0.27 

0.24 

0.40 

0.28 

0.26 

0.26 

0.28 

0.45 

0.29 

0.29 

0.26 

0.27 

0.50 

0.29 

0.27 

0.25 

0.29 

0.55 

0.30 

0.26 

0.27 

0.28 

0.60 

0.30 

0.28 

0.29 

0.30 

0.65 

0.32 

0.28 

0.27 

0.30 

0.70 

0.32 

0.28 

0.29 

0.29 

0.75 

0.31 

0.29 

0.30 

0.29 

0.80 

0.27 

0.28 

0.28 

0.27 

0.85 

0.29 

0.27 

0.26 

0.26 

0.90 

0.27 

0.28 

0.28 

0.25 

0.95 

0.26 

0.26 

0.26 

0.27 

1.00 

0.28 

0.21 

0.26 

0.25 
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APPENDIX  G 

ELECTROCHEMICAL  MEASUREMENTS 

The  physical  and  chemical  properties  of  the  elctrochemical  solution 
are  given  in  this  appendix. 

Physical  Properties 

All  the  electrochemical  measurements  required  the  determination  of 
the  kinematic  viscosity,  density,  concentration  of  active  specie  and 
mass  diffusivity.  The  viscosity  was  measured  after  each  run  with  a 
No.  50  Otswald  viscometer.  The  density  was  determined  by  hydrometer 
immersion.  The  concentration  was  obtained  using  standard  titration 
for  potassium  iodide.  The  diffusion  coefficient  for  the  iodine  was 
calculated  from  the  following  correlation  (see  Zilker  [77]), 

log  D  =  -1.07291  log1(J  v  -  7.15278  (G.l) 

where  v  is  the  kinematic  viscosity.  Table  G.l  summarizes  the  physical 
and  chemical  properties  of  the  electrolyte. 


228 


Table  G.l. 

Physical  and  Chemical  Properties  of  Electrolyte 


Property 

Viscosity  (P) 

0.00858 

3 

Density  (g/cm  ) 

1.023 

2 

Kinematic  Viscosity  (cm  /s) 

0.00839 

2 

Hass  Diffusivity  (cm  /sec) 

1.188  x  10-5 

Schmidt  Number 

707 

KI  Concentration  (moles/liter) 

0.103  -  .2 

*2  ^3^  Concentration  (moles/liter) 

0.00131  -  0.00146 

Temperature  (C°) 

26 
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APPENDIX  H 
ANALYSIS  OF  RAW  DATA 


In  this  appendix  the  method  used  in  analyzing  the  raw  data  is 
outlined.  All  electrodes,  are  assumed  to  have  the  same  area.  The 
instantaneous  current  measured  at  the  jth  electrode  is 


E. 


where  E  *  (V  -  V  N .  Here  V  is  the  instantaneous  measured  voltage 
3  °  aPP)  ° 

and  V  is  the  applied  voltage.  The  term  is  the  feedback  resistance, 
app  f 

The  average  current  at  the  jth  electrode  is  given  by 


I . 
J 


n 

l 


i=l 


R^n 


(H.  2) 


where  E^  is  the  ith  voltage  sampled  by  the  A/D  converter  and  n  is  the 
total  number  of  samples  gathered.  The  average  shear  stress  at  the 
jth  electrode  is  obtained  from  equation  (4.5) 


(t  )  -  (T 

wd 


V  E1  .3 

l  j 

i=l 


nR, 


(H.  3) 


3  (E1>' 
=  c  — -V 


(H.  3) 


r: 


n 

l 

i=l 


Ej/n. 


where  E . 

J 
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The  wave  averaged  shear  stress  is  obtained  by  summing  (t  ) .  over 


all  j, 


<T 


m 

£  (t  )  . 
•  T  w  J  1 
>  =  1=1  d 


w , 


m 


(H.  5) 


where  m  is  number  of  electrodes.  The  dimensionless  quantity, 

(t  /  x  )  =  (t  )./<x  >,  is  then 

w-j  w  wd  1  wd 


T 
W 


(E.)3 

- J _ 


t-  m  o 

l  (E.)  /m 


tj=l 


(H.  6) 


(EJ3 


(E3) 


w 


(H.  7) 


where 


I  <v3 

- 


m 


(H.  8) 


The  local  intensity  of  the  fluctuating  shear  stress  is  found  using 
equation  (4.6)  to  be 


rr  i  (E .  -  E . ) 


T  ,  . 

W;] 


=  3  _J _ 1 


(H.  9) 


231 


The  local  root  mean  square  value  of  the  fluctuating  shear  stress  is 


j 


-  3 


E. 

3 


(H. 10) 


and  the  absolute  root  mean  square  value  of  the  fluctuating  shear  stress 
is  obtained  by  multiplying  equation  (H.9)  by  equation  (H.6) 


2  —2  1/2  _o 

(ET  -  ET)  E. 

=  3  _ 3_  _ 3— 


(V 


<s\ 


(H. 11) 


X 

S 

X 

where  s1  =  3uf/3y  and  S  *  3U/3y  are  respectively,  the  dimensionless 
x  x 

fluctuating  and  wave  averaged  velocity  gradient  evaluated  at  the  wall. 
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APPENDIX  I 


LEAST  SQUARES  ANALYSIS 


The  surface  shear  stress  data  are  fitted  with  a  Fourier  series. 

The  computer  program  Data  is  set  up  to  calculate  the  first  four  harmonics 
if  necessary,  however,  onlyone  harmonic  was  used  in  fitting  the  data.  The 
Fourier  coefficients  are  obtained  by  performing  a  least  squares  analysis 
as  follows: 

let  f(x)  be  a  function  that  has  a  period  2Tt/a.  Then  the  Fourier  series 
representation  of  f (x)  is 


_  N 

f(x  )  *  f  +  £  (a  cos  (anx  )  +  b  sin  (anx  ) ) 

J  i  ^  j  ^  J 


(1.1) 


J 


l  fOO 
1  J 


(1.2) 


J 


where  x^  is  the  location  of  jth  electrode.  The  spacing  between  the 
electrodes  is  assumed  constant  in  the  analysis  that  follows. 

The  error  function  E  is  defined  as 

J 


(1.3) 


where 


2 


(1-4) 
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The  constants  a  and  b  are  found  by  minimizing  the  error  function  with 
^  n  n 

respect  to  a  and  b 
r  n  n 


3E 

9a 

n 


0, 


3E 

9b 


=0  n=l,  2, 


(1.5) 


Substituting  equation  (1.4)  into  (1.3)  and  using  equation  (1.5) 
yields 


Equations  (1.6)  and  (1.7)  represent  2N  equations  in  2N  unknowns, 

a  and  b  (n  =  1,  . .  . ,  N) .  The  function  f (x) ,  can  be  written  as  follows 
n  n 


_  N  o  o  1/  2 

f(x)  =  f  +  y  (a  +  b  )  cos  (anx  +  0  ) 
“  n  n  n 


where 


0  =  tan  1  (-b  /a  ) 

n  n  n 


(1.8) 


(1.9) 
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APPENDIX  J 

CORRECTION  OF  INTENSITY  DATA 


In  this  appendix  the  measurements  of  /sTZ  / S  are  corrected 

v  x  x 

for  both  frequency  and  spatial  averaging  as  suggested  by  Mitchell  and 

Hanratty  [47].  The  results  are  summarized  in  Tables  J.l  and  J.2, 

where  L+  is  the  effective  electrode  length,  W  =  2 II  (0.009)  is  the 

m 

median  frequency  of  in  pipe  turbulence.  Sc  is  the  Schmidt  number 
and  A  is  the  circumferential  scale  for  mass  transfer  fluctuations. 
Mitchell  and  Hanratty  have  suggested  that  for  pipe  turbulence  A  *  12. 


They  have  presented  curves  of  1/A2  vs.  (L+W+Sc1//2)  and  /  k ' 2  / 


m 


accounts 


2  1/2 

vs.  (1/A)  where  (1/A  )  is  the  correction  factor  which 
for  frequency  averaging  and  / k'2  /  /k^2  is  the  correction  factor  which 
accounts  for  spatial  averaging.  Here  k'2  is  the  corrected  intensity 
of  the  mass  transfer  fluctuation  and  k/  is  the  measured  intensity  of 
the  mass  transfer  fluctuations. 


It  should  be  noted  that  the  values  of  /st2/S  listed  in  this 

v  x  x 

appendix  refer  to  the  wave  induced  components  of  the  fluctuating 


shear  stress. 


Table  J.l 


Correction  of  Fluctuating  Shear  Stress  Measurements 


a 

L+ 

+  +3/2  1/2 

L  W  Sc 

m 

L/A 

A’2 

1 

P 

v  m 

7 

0.00787 

6.442 

3.355 

.5369 

- 

1.09 

0j 00621 

8.162 

4.25 

.6802 

- 

1.22 

0.00565 

8.974 

4.674 

.748 

- 

1.26 

0.0045 

11.267 

5.868 

.939 

1.07 

1.29 

0.00367 

13.816 

7.915 

1.151 

1.11 

1.40 

0.00297 

17.07 

8.9 

1.422 

1.15 

1.48 

0.00287 

17.66 

9.201 

1.472 

1.155 

1.5 
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Table  J.2 

Correction  of  Fluctuating  Shear  Stress  Measurements 


a 

/Sx 

/<  M 

/si2  i 

1 

sx 

\  /W 

s*  A 

X 

V  m 

A 

1.00787 

0.0462 

0.0462 

0.0483 

0.048 

1.00621 

0.0482 

0.0482 

0.053 

0.053 

1.00565 

0.0465 

0.0465 

0.0528 

0.0528 

1.00435 

0.0366 

0.0392 

0.0416 

0.0445 

1.00367 

0.0305 

0.0338 

0.0374 

0.0397 

1.00297 

0.0211 

0.0243 

0.025 

0.0295 

1.00287 

0.0230 

0.0265 

0.0280 

0.0325 
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APPENDIX  K 

NUMERICAL  FORMULATION  OF  WAVE  INDUCED  EQUATIONS 


In  this  appendix  equations (3. 13)- (3 . 17)  are  set  up  in  a  suitable 
form  for  numerical  solution  using  the  techniques  described  in  Chapter  5. 

The  differential  equation  describing  the  wave  induced  flow  using 
the  zero  equation  turbulence  models  is  a  fourth  order  equation.  Since 
the  dependent  variable,  the  stream  function,  F,  is  a  complex  quantity 
the  differential  equation  is  solved  as  an  eighth  order  system.  This 
is  achieved  by  splitting  F  into  its  real  and  imaginary  parts  and 
introducing  the  following  variables. 


=  F. 


R 


^  =  F, 


F 1 
R 


F  * 
I 


fr 


FI 


♦  a  -  F 


R 


T7n  * 

FI 


(K.  1) 


Equation  (3.13)  can  be  written  as  follows: 


F,m  =  aQ(y)F  +  a^F'  +  a2(y)F"  +  a3(y)F,,T  +  g(y)  (K.2) 


where  a^Cy),  a^(y),  a2(y)  and  a^(y)  are  the  complex  coefficients  of 
the  derivatives  of  the  stream  function  and  g(y)  is  the  complex  non- 
homogeneous  term.  Equation  (K.2),  can  be  expressed  in  real  and 
imaginary  form  as  follows. 
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F""  =  f""  +  iF"" 
R  I 


(K.3.1) 


a0RFR  aOrFI  +  31RFR  al].FI  +  a2RFR  a2IFI  +  33rFr'  “ 


Q  J*  |  ^  p  _|_  Q  pi  .p.  £  pi  ^  pH  ^  pM  _j_  ^  pill  pill 

L  °I  R  °R  1  FI  R  XR  1  2I  R  2R  1  3I  R  3R  1 


+  gR  +  igT 


(K. 3 . 2) 


where 


a  = 


= 


= 


a  tia 
R  I 


a  +  i  a 
R  I 


a0  ■' 


r  a„  +  i 
JR 


g  = 


gR  +  i8I 


(K.4.1) 

(K.4.2) 

(K.4.3) 

(K.4.4) 

(K.4.5) 


Equation  (K. 3.2)  can  be  written  in  a  form  similar  to  equation  (5.61), 


_d_ 

dy 


■*l‘ 
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The  boundary  conditions  are 

F  =  0  F’  =  0  at  y  =  0 

F  =  U  F !  =  U 1  at  large  y. 


In  matrix  notation  these  conditions  are:  at  y  =  0, 
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Since  the  number  of  boundary  conditions  at  the  wave  surface  and 
at  large  y  are  both  equal  to  two,  there  is  no  preferred  direction  of 
integration.  In  the  computations  carried  out  in  this  thesis,  the 
direction  of  integration  is  toward  the  wave  surface.  Hence  y  =  0  =  b 
and  y  — >  y  =  a. 

J  J  oo 


The  particular  solution  is  generated  by  solving  the  system  of 
equations  (K.5)  subject  to  the  boundary  condition 


U 

U’ 

0 

0 

0 

0 

0 

0 


at  large  y. 


(K.9) 


The  required  solution  to  equation  (K.3.2)  could  be  obtained  by 
generating  a  particular  solution  satisfying  equation  (K.9)  and  eight 
independent  homogeneous  solutions.  The  required  solution  would  then  be 
obtained  by  forming  a  linear  combination  of  the  eight  homogeneous 
solutions  and  the  particular  solution.  However  by  suitable  choice  of 
the  homogeneous  solutions  at  y^  it  is  necessary  to  generate  only  four 
homogeneous  solutions.  This  is  demonstrated  below. 
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Let  the  eight  homogeneous  solutions  at  satisfy 
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The  required  solution  is, 


F  =  F  +  c-  Fu  +  c  F  +  c  F  +  c  F  +  c  F__  4*  c  F 

p  1  2  H2  3  H3  4  5  H5  6  H6 


-f  c  F  4*  c  F 

7  H7  8  H8  ‘ 


(K. 11) 


The  above  combination  must  satisfy  the  boundary  conditions  at 
yw,  equation  (K.8),  hence  using  these  conditions, 

"00001000  '' 

00000100 
00000010 
00000001 


In  order  for  the  above  equations  to  be  satisfied  c,- ,  c^,  c^ 
and  Cg  must  be  zero.  Therefore,  it  is  necessary  to  generate  only 
four  independent  homogeneous  solutions.  The  constants, c^  -  c^,are 
determined  once  Fu  -  ¥  and  F  have  been  integrated  to  y  =  0,  by 

H1  \  P 

satisfying  the  boundary  conditions  at  y  =  0. 

The  K-e  Model  equations  (3.72)— (3.77.2)  was  solved  using  the  same 
technique  outlined  above.  The  K-e  Model  requires  the  solution  of  an 
eighth  order  complex  differential  equation  which  is  then  converted 
into  a  sixteenth  order  real  system.  This  is  achieved  by  splitting 
the  stream  function,  F,  the  kinetic  energy,  k,  and  the  dissipation 
rate,  e,  into  real  and  imaginary  parts  and  introducing  the  following 


variables: 
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(K. 13) 


The  K-e  Model  equations  (3.13),  (3.72)  and  (3.73)  can  be  organized 
in  a  similar  manner  to  equation  (K-3) 


F""  =  aQ(y)F  +  a1(y)F*  +  a2(y)F"  +  a^'" 


+  a^(y)k  +  a5(y)k'  +  ag(y)e  +  a7(y)e'  +  g  , 


(K. 14) 


k"  =  bQ(y)F  +  b1(y)F'  +  b4(y)k  +  b5(y)k 


f  b^ (y) e  +  b7(y)e'  +  f  , 


(K. 15) 


e"  -  cQ(y)F  +  c^(y)F*  +  c4(y)k  +  c^(y)k* 


+  c&(y)e  +  c7(y)e?  +  h  , 


(K. 16) 


where  the  a^(y),  b^(y)  and  c^(y)  are  the  complex  coefficients  such 


that 
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(K.17.1) 


(K.17.2) 


(K.17.3) 


The  system  of  equations  (K. 14)-(K-17)  can  be  written  in  a  form 
similar  to  equation  (5.6) 
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where 


(K. 19) 
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The  boundary  conditions  at  y  =  0  cannot  be  applied  directly  since 
equations  (K.15)  and  (K.16)  are  singular.  Therefore,  this  boundary 
condition  is  applied  at  a  small  distance,  y  =  <5^,  from  the  wave  surface. 
In  order  to  apply  the  boundary  condition  a  Taylor  series  approximation 
is  made : 
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where 


(K. 24) 


(K. 25) 


(K. 26) 


(K.27) 
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The  required  solution  to  equations  (K. 18)- (K. 28)  is  obtained  by 
generating  a  particular  solution  and  eight  homogeneous  solutions.  The 
particular  solution  is  chosen  to  satisfy  the  known  boundary  conditions 
at  y  .  The  initial  vectors  needed  to  generate  the  eight  homogeneous 
solutions  are  obtained  in  a  similar  manner  described  for  equation  (K.5). 
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APPENDIX  L 

RELATIONSHIP  BETWEEN  TURBULENT  ENERGY 
AND  VELOCITY  GRADIENT  AT  THE  WAVE  SURFACE 

In  this  appendix  the  wave  induced  kinetic  energy  is  related  to 
the  root  mean  square  level  of  the  fluctuating  velocity  gradient, 


All  quantities  are  made  dimensionless  with  wall  parameters. 


The  fluctuating  velocity  close  to  the  wave  surface  is  related  to  the 
fluctuating  velocity  gradient  as 
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Now  R^  can  be  decomposed  into  a  wave  averaged  and  fluctuating 
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Hence: 


xx 


(L.5) 
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a  r 


where 


2 (-R  )1/2 

XX 


xx  lax  . 

e  is  the  wave 


induced  component  of,  y  J .  The  turbulent  kinetic  energy 
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K  is  related  to  u’  as  follows 
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Since  -  < -  and  -  << — -  in  the  neighborhood  of  the  wall 
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(L.7) 


The  turbulent  kinetic  energy  can  be  written  in  terms  of  a  wave  averaged 
and  fluctuating  component  as 
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Equations  (L.5)  and  (L.9)  gives  the  desired  result  . 
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where  0^  is  the  phase  angle  of  the  wave  induced  variation  of 
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NOMENCLATURE 


a 

ad 

A 

A 

A 


m 


I 

A 

k 


Dimensionless  wave  amplitude 
Wave  amplitude  (in.) 

Van  Driest  parameter,  see  equation (3. 47) 

Average  value  of  Van  Driest  parameter 

Complex  function  related  to  the  wave  induced  component 
of  the  Van  Driest  Parameter,  see  equation  (3.48.1) 

2 

Electrode  area  (cm  ) 

Constant  associated  with  K-e  Model  defined  in  equation  (3.62) 

Constant  associated  with  K-e  Model  defined  in  equation  (3.62) 

Constant  associated  with  K-e  Model  defined  in  equation  (3.64) 

Constant  associated  with  K-e  Model  defined  in  equation  (3.63.2) 

Constant  associated  with  K-e  Model  defined  in  equation  (3.66) 

Concentration  of  diffusing  species  in  bulk  (moles/liter) 
Diffusion  coefficient  for  mass  (cm/sec) 

Damping  function,  see  equation  (3.35) 

Function  associated  with  K-e  Model 

Function  associated  with  K-e  Model 

Dimensionless  complex:  function  related  to  the  wave  induced 
stream  function 

Faradays  constant  (coulombs/equivalent) 

Ch&nnal  height  (in.) 

Dimensionless  metric  in  the  x  direction 
Dimensionless  metric  in  the  y  direction 
Cathode  current  (amps) 

Dimensionless  complex  function  related  to  the  wave  induced 
turbulent  kinetic  energy 
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kft 


k 

c 

K 


K 


K 

m 

£ 


a 

£ 


£ 

o 


L 

e 

n 

e 


P 


p' 

A 

P 


Coefficients  defined  in  equation  (3.47) 

Lag  parameter  defined  in  equation  (3.53) 

Lag  parameter  defined  in  equation  (3.56) 

Lag  parameter  defined  in  equation  (3.82) 

Dimensionless  turbulent  kinetic  energy 

Dimensionless  wave  averaged  turbulent  kinetic  energy 

Mass  transfer  coefficeint  defined  in  equation  (4.1) 

Dimensionless  mixing  length;  modified  for  streamline 
curvature 

Dimensionless  complex  function  related  to  the  wave 
induced  mixing  length;  modified  for  curvature 

Dimensionless  plane  shear  mixing  length 

Dimensionless  complex  function  related  to  the 
wave  induced  plane  shear  mixing  length 

Equivalent  length  of  electrode  (cm.) 

Number  of  electrons  transferred  in  reaction 

Dimensionless  instantaneous  pressure 

Dimensionless  fluctuating  pressure 

Dimensionless  complex  function  related  to  the 
wave  induced  pressure 


P  Dimensionless  time  averaged  pressure 

—  Dimensionless  value  of  pressure  gradient 

X  in  x  direction 


Dimensionless  average  value  of  pressure  gradient 
in  x  direction 


3P 

3y 


Dimensionless  value  of  pressure  gradient 
in  y  direction 


3P 

ay 


Dimensionless  average  value  of  pressure  gradient 
in  y  direction 


2 

q 


Dimensionless  trace  of  Reynolds  stress  tensor 


r 


ij 


Dimensionless  function  related  to  the  wave  induced 
components  of  the  turbulent  stress  tensor 


Cfl  > 
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R 

c 

Re 


R 


s 

X 


U 

U 


Dimensionless  components  of  the  turbulent  stress  tensor 

Dimensionless  average  components  of  the  turbulent 
stress  tensor 

Dimensionless  radius  of  curvature  of  streamlines 
Channel  Reynolds  number  based  on  and  channel  half  width 
Feedback  resistance  (ohms.) 

Dimensionless  curvature  Richardson  number  defined 
on  equation  (3.80) 

Dimensionless  complex  function  associated  with  the 
wave  induced  Richardson  number 

Dimensionless  composite  term  representing  turbulent 
stress  terms,  see  equation  (3.14) 

Dimensionless  fluctuating  component  of  velocity 
gradient  at  the  wall 

Dimensionless  complex  function  related  to  the  wave 
induced  rate  of  stress  tensor 

Dimensionless  components  of  the  rate  of  strain  tensor 

Dimensionless  average  rate  of  strain 

Dimensionless  velocity  gradient  at  the  wall 

Dimensionless  wave  averaged  velocity  gradient  at  the  wall 

Dimensionless  instantaneous  velocity  in  x  direction 

Dimensionless  fluctuating  component  of  turbulent 
velocity  in  x  direction 

Dimensionless  complex  function  related  to  the  wave 
induced  velocity  in  the  x  direction 

Friction  velocity  (ft/sec) 

Dimensionless  time  averaged  velocity  in  the  x  direction 
Dimensionless  wave  averaged  velocity  in  the  x  direction 
Bulk  average  channel  velocity  (cm/sec) 

Dimensionless  instantaneous  velocity  in  the  y  direction 


V 
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vr  Dimensionless  fluctuating  component  of  turbulent 

velocity  in  the  y  direction 

v  Dimensionless  complex  function  related,  to  the  wave 

induced  velocity  in  the  y  direction 

V  Dimensionless  time  averaged  velocity  in  the  y  direction 

V^pp  Applied  voltage  (volts) 

Vq  Instantaneous  measured  voltage  (volts) 

x  Dimensionless  boundary  layer  coordinate  along  wave 

y  Dimensionless  boundary  layer  coordinate  perpendicular 

to  wave 

Greek 

a  Dimensionless  wave  number 

Coefficient  defined  in  equation  (3.80) 

Y  Coefficient  defined  in  equation  (3.91) 

6  Boundary  layer  thickness 

6  ,  Small  distance  from  wave  surface  defined 

in  equation  (3.78.1) 

/s 

e  Dimensionless  complex  function  associated  with  wave 

induced  turbulent  dissipation  rate 

e  Dimensionless  turbulent  dissipation  rate 

e  Dimensionless  average  turbulent  dissipation  rate 

0  Phase  angle  (°) 

K  Curvature  of  wave  surface 

k  Von  Karman  constant  used  in  equations  (3.35), 

(3.41),  (3.43),  (3.50),  (3.58),  (3.71),  (3.86) 

and  (3.69) 

X  Wave  length  (in*) 

p  Viscosity  of  fluid  (poise) 

2 

v  Kinematic  viscosity  (cm  /sec) 


H  >  H  ~0  <L\  <1> 
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v 


t 

t 


t 


2 

Turbulent  viscosity  (cm  /sec) 

Complex  function  related  to  the  wave  induced  turbulent 
viscosity  (cm/ sec) 

Average  turbulent  viscosity  (cm/sec) 

Coefficient  defined  in  equation  (3.89) 

3  ‘  * 

Density  of  fluid  (gm/cm  ) 

Dimensionless  shear  stress 

Dimensionless  complex  function  related  to  the  wave 
induced  shear  stress 

Dimensionless  wave  averaged  shear  stress 

Dimensionless  wall  shear  stress 

Dimensionless  wave  averaged  wall  shear  stress 

Dimensionless  fluctuating  shear  stress  in  x  direction 

2 

Instantaneous  wall  shear  stress  (dynes/cm  ) 

2 

Time  averaged  wall  shear  stress  (dynes/cm  ) 

Fluctuating  component  of  turbulent  shear  stress  in 
the  x  direction  (dynes/cm^) 

Dimensionless  stream  function 
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